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ABSTRACT 


This research used Huygens-Fresnel wave optics computer simulations to 
investigate the effects of high turbulence strength and inner scale on the 
normalized irradiance variance and coherence length of electromagnetic waves 
propagating through a turbulent atmosphere. These investigations developed 
several guidelines for validity of propagation simulations employing a numerical, 
split-step, Huygens-Fresnel, method, and within these guidelines, considered five 
types of turbulence spectrum inner scale: (1) zero inner scale, (2) Gaussian inner 
scale, (3) Hill's and (4) Frehlich's viscous-convective enhancement inner scales, 
and (5) turbulence spectrum truncation from the discrete grid representation. The 
simulation results showed that the normalized irradiance variance generally 
decreased (-30%) below the zero inner scale values in the Rytov regime with 
increasing inner scale size, but increased monotonically in the saturation regime, 
and agreed within 2% of the Rytov-Tatarski predictions at low turbulence strengths. 
The E-field coherence length in a spatially confined beam, with either spherical or 
plane wave divergence and zero inner scale, followed the Rytov-Tatarski-Fried 
predictions in the Rytov regime, but departed from the theory in the saturation 
regime. Increasing inner scale size modified this finite beam behavior by raising 
the coherence length (up to -50%) in the saturation regime. 
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I. INTRODUCTION 


Turbulence in a stratified fluid causes random inhomogeneities in 
temperature and the index of refraction that scatter and diffract a wave 
propagating through such a medium. Analytical perturbation techniques 
developed over the past three decades cannot account for the variations 
observed experimentally in the amplitude and phase distributions of a 
propagating wave when turbulence levels are high and/or propagation paths are 
long. Over the past fifteen years, numerical wave optics codes based on the 
Huygens-Fresnel principle have been developed to address these situations. 
This research extended these codes to include inner scale in the spectrum of 
refractive index fluctuations and to examine the coherence length and the 
effects of inner scale in the Rytov regime (low turbulence strengths and/or short 
propagation paths) and in the saturation regime (high turbulence strengths 
and/or long propagation paths). 

As a first step, extensive analysis and testing developed guidelines for 


validity of computer simulations employing Huygens-Fresnel propagation over 
multiple steps (split-step method). These guidelines include: 



For an NxN grid and propagation distance L, the grid element size should 
be Ax = yTTTlV . 

The maximum turbulence strengths C„ 2 and propagation distances L for 
spherical wave propagation with an NxN grid should satisfy 
Po mai S 0.1 N 09 , where po = 0.497 C„ J k 7 ' 6 L” ,e . 


• Equivalently, the E-field coherence length r 0 should satisfy 

r 0 a 2.5 Ax, where r 0 represents Fried's coherence length (for 
spherical waves). 

• Use > 30 phase screens/steps for each propagation. 

• Use > 30 propagation realizations to get representative statistics. 

• Phase screens require low spatial frequency correction to gain 5% 
accuracy in normalized irradiance variance and as much as 30% in 
coherence length. 

• Half width at half maximum of the atmospheric MTF and an iterative fit r 0 
provide the most stable parameterizations of the E-field coherence length. 

• Telltale signs of aliasing include a fine-grained irradiance pattern, a boxed 
perimeter to the irradiance pattern, and peaking of the energy toward the 
center of the computation grid. 

These investigations also examined four choices of E-field source function and 
refined the methods of normalized irradiance variance and E-field coherence 
length calculation. 

These simulations incorporated five types of turbulence spectrum inner 
scale: (1) zero inner scale, (2) Gaussian inner scale, (3) Hill's and (4) Frehlich’s 
viscous-convective enhancement inner scales, and (5) turbulence spectrum 
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truncation from the discrete grid representation. For the more physically 
realistic viscous-convective enhancement inner scale, the computer simulations 
provided the following results: 

• The Hill and Frehlich parameterizations performed almost identically, 
giving less than 3% difference in normalized irradiance variance over the 
Rytov and saturation regimes. 

• In the Rytov regime, the normalized irradiance variance for an 
approximately spherically diverging beam increased (-30%) and then 
decreased (-30%) compared to the zero inner scale values as the inner 
scale size increased, and the E-field coherence length r 0 decreased 
slightly (-5%) compared to the zero inner scale coherence length. 

• In the saturation regime, increasing inner scale size gave monotonically 
increasing normalized irradiance variance and increased the coherence 
length (up to -50%) compared to the zero inner scale case. 

The effect of the Gaussian inner scale on normalized irradiance variance and 

coherence length was also investigated. 

These investigations examined the behavior of the E-field coherence 

length for E-fields (beams) that were constrained in lateral extent to the grid 

size but whose divergence approximated that of spherical waves, plane waves, 

and intermediate beam waves. The computer simulations showed that: 

• For the Rytov regime, the spherical and plane wave cases gave 
coherence lengths within 5% of the Rytov-Tatarski-Fried predictions. 

• In the saturation regime, the E-field coherence lengths for the spherical 
wave approximation decreased -25% below the theory, while the 
coherence lengths for the plane wave approximation varied within 
-5%/+15% of the theory. 
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• The E-field coherence lengths for the beam wave cases were spaced 
between the spherical and plane wave values for the Rytov regime, but 
increased toward the spherical wave values in the saturation regime. 

The spatially confined beams used in these simulations showed a departure in 

behavior of the E-field coherence length in the saturation regime from the 

predictions of first order perturbation theory. 

The organization of this dissertation generally follows the preceding 

summary of major points. Chapter II summarizes the theory of wave optics and 

the implementation of the Huygens-Fresnel propagation in a computer 

simulation. It then provides an overview of the Rytov-Tatarski theory for the 

effect of atmospheric turbulence, including inner scale, on the normalized 

irradiance variance and E-field coherence length as applied to spherical wave 

propagation. Chapter ill discusses the detailed choices of physical and 

simulation parameters and develops the guidelines to ensure validity of the 

simulation. Chapter IV presents the computer simulation results, first for the 

inner scale effects on normalized irradiance variance, then for the behavior of 

the E-field coherence length in the Rytov and saturation regimes, and for the 

effect of inner scale on coherence length. 
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II. BACKGROUND 


A. CHAPTER OVERVIEW 

This chapter summarizes first the theory of electromagnetic wave 
propagation and the Huygens-Fresnel solution to the scalar wave equation and 
recasts the Huygens-Fresnel solution into a form utilizing fast Fourier 
transforms easily implemented in a computer simulation. It then summarizes 
the Rytov-Tatarski solution to the scalar wave equation which used first order 
perturbation theory and statistical techniques to describe the spatial variations 
of the irradiance of an electromagnetic wave propagated through a turbulent 
medium. The Kolmogorov spectrum of refractive index fluctuations is 
introduced as well as five types of inner scale which modify the high spatial 
frequency portion of this spectrum. The chapter concludes by describing 
Fried's method for parameterizing the coherence length of the E-field for a 
spherically-diverging wave propagating through a turbulent medium. 

B. WAVE PROPAGATION 

Maxwell's equations describe of the propagation of electromagnetic 


waves through a turbulent medium. Assuming a locally homogeneous, 
isotropic, and linear medium, 



V ■ cE 


(1) 


= P- 
V ■ B = 0, 


V 


. as 


( 2 ) 

(3) 


V x H = 

dt 


(4) 


where hatted ( A ) quantities represent vectors. At the frequencies and field 
strengths of interest, the medium may be assumed to have zero local free 
charge density p and zero (or negligible) conductivity o. Laser beams in the 
atmosphere satisfy these conditions. Expanding the left-hand side of Eq. (1) 
and combining the curl of Eq. (3) and the time derivative of Eq. (4) gives a 
vector wave equation for the E-field 


V 2 £ + V(£ • V In e) - n c = 0. (5) 

The middle term represents depolarization of the E field and is negligibly small 
for propagation through a turbulent atmosphere (Tatarski, 1961; Lawrence and 
Strohbehn, 1970). Neglecting the middle term simplifies the vector wave 
equation to 


P 


PE 

dt 2 


V 2 £ - 


( 6 ) 



The Fresnel-Kirchoff diffraction theory provides an approximate scalar 
solution to this vector wave equation. Following Hecht (1987), the spherical 
wave solution is 

E(U) = ~ e m ' ut) , ( 7 ) 

where E 0 is a constant and k = 2k JX . Using Eq. (7) with the time dependence 
factored out, the Kirchoff integral theorem becomes 

a ■ i$ls £ r v£ «>- < * s - lls ^1-■“*] • «•» 


which must be evaluated over a surface S enclosing the field point P. Figure 
(1) illustrates the relationship between the distances £ and ry If the wavelength 
X « £,, q , Eq. (8) reduces to the Fresnel-Kirchoff diffraction formula 


Ep 


Eq i e W'n) 
X JJs 5 n 


K(Q) dS, 


(9) 


where K(0) represents the obliquity factor. 

Integrating Eq. (9) over the half sphere S (shown in Fig. (1)) as the 
radius approaches infinity, the electric field amplitude is appreciable only over a 


finite aperture area A of the hemisphere's flat side. Then Eq. (9) reduces to the 
more familiar Huygens-Fresnel expression 




Figure 1 Fresnel-Kirchoff diffraction geometry showing relationship between 
source, aperture A, surface S, and field point P. 


£p = j / l A E a /C(0) dA. (10) 


The complex aperture function E A contains the source spherical wave factor, 
exp(ik£)/£, . The E-field at a point P is a function of the E-field over a finite 
aperture A. 
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The Huygens-Fresnel formulation , Eq. (10), serves as the basis for 
calculating the propagation of a wave through turbulence. This formula can 
now be cast into a form easily implemented in a propagation computer 
simulation. Applying the Cartesian coordinate system of Fig. (2), the Huygens- 
Fresnel principle becomes, 

-i rr +■ \f - p| z 

Efrz) = ff E(p,0) e m dA. (11) 

A Jz 2 + \f - p | 2 



Aperture 


Figure 2 Propagation coordinate system showing relations between the 
aperture variables and the field variables. 



t represents the position vector of the observation point P in the x-y plane and 
p represents the position vector of the radiating point in the aperture plane. 

The paraxial approximation assumes that 

|?| < z and |p| « z, 02) 

where z measures the longitudinal distance of the point P from the aperture. 
Consequently, the obliquity factor K(0)=(1+cos0)/2 becomes approximately 1 
and the denominator of the integrand is approximately z (Goodman, 1968), 

E{t,z) E( p,0) I'-pI* dA . (13) 


Following Roberts (1986), the Fresnel approximation assumes that the lateral 
displacement between the radiating point and the observation point 
| t - p | « z. Then, using the lowest order terms for the series expansion 
of the complex exponential's argument, Eq. (13) becomes 


£f(Az) ~-^\d P E(p,0) e 


iif<fp E(P.O) * 




(14) 


- e k J dp E{ p,0) e x 
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The factor exp(-27tiz/X,) represents the change in phase of the optical wave from 
the center of the aperture plane to the plane of the observation point P. It 
applies to the whole E-field and depends only on the propagation distance z. It 
does not affect the phase variation across the E-field nor the amplitude of the 
E-field, and will be dropped in the following expressions for simplicity. 

Equation (14) can be recast into a form utilizing Fourier transforms. 
Expanding the aperture field E(p,0) in Eq. (14) using the Fourier transform 
identity 

£(p,0) = [df e 2n/, p [dp' e-*" 1 *' E(p',0), (15) 

Eq. (14) becomes 

E(f,zy --L [dp [ [dfe 2 ^ fdp' e 2 ""*' E( p',0) ] (16) 

(Note: these equations use spatial frequency / (m' 1 ) instead of spatial 
wavenumber i = 2nf (rad/m) because the discrete Fourier transform 
computer subroutines used in these investigations are written with spatial 
frequency f.) Making the change of variables f e t - p gives 


11 



( 17 ) 


E[f,z) = --L f(-a') [ e > 2 nHr-r') fffri e-^h'E(p',Q) | e** 1 ' 1 , 

= j- jdf o a « ft fdp' E(p',0) [ jdf‘ e~ an,f ' ]. 

The last integration over 6f is the Fourier transform of a Gaussian function, 

I 0 , e - a «rr . iXz (18) 

Substituting p for p', the E-field of Eq. (17) becomes 

E(t,z) = fd?e a * ft fdp E(p,0) e (19) 

Equation (19) may be symbolically written as 

E(f,z) = IFT[ g-kizft* Fr[E(p,0)]], (20 > 

where FT and IFT represent the two-dimensional Fourier transform and inverse 
Fourier transform, respectively. Equation (20) expresses the E-field at a 
propagation distance z in terms of Fourier transforms of the E-field at z = 0 
represented in a Cartesian coordinate system. This form of the Huygens- 
Fresnel principle is useful for a simulation that subdivides the propagation path 
into a sequence of short Fresnel propagation steps. 
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C. EFFECTS OF TURBULENCE ON THE PROPAGATED WAVE 
1. Irradiance Variance and Inner Scale 

The next step in modeling the propagation of an E-field through a 
turbulent medium requires a statistical characterization of the turbulence and its 
effect on the spatial statistics of the propagating wave. To begin understanding 
the effect of turbulence upon the wave, consider a spherical wave incident upon 
a medium with randomly varying index of refraction, as shown in Fig. (3). For 
small variations in the index of refraction about the mean value, scattering 
occurs predominantly in the forward direction (Tatarski, 1961). The regions of 
index of refraction fluctuation accelerate/retard portions of the wavefront over 
short propagation distances and cause variations in phase across the reference 
spherical wavefront. Subsequent diffraction and interference then create 
variations in irradiance across the wavefront. Compared to propagation through 
zero turbulence (Fig. (4)), a spherically diverging beam now exhibits significant 
variations in phase and amplitude (irradiance) (Fig. (5)). 

A statistical description of the refractive index fluctuations is 
needed to calculate the E-field variations. The index of refraction, n, in air 
depends upon the temperature T and pressure P (Tatarski, 1961). At visible 
wavelengths, 
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Figure 3 Wave propagating through a medium having random inhomogeneities 
in index of refraction. 


n - 1 . 79 xio* PJOWgrtL . (21) 

T (a) 

The temperature profile in the atmosphere possesses significant stratification, 
as shown in Fig. (6) (Walters, 1994). Velocity shear between different layers in 
the atmosphere causes turbulence which disrupts the interface between these 
layers, mixing regions of different temperatures. The velocity fluctuations 
generated by the turbulence at the interface between layers generally follow the 
Kolmogorov spectrum, 


14 










Figure 4 Irradiance plot for a spherically diverging beam propagated through 
zero turbulence. 
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Figure 5 Irradiance plot for spherically diverging beam propagated through 
turbulence. 
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rfcj 

nurarma^M 
0*» FWi - Sa#NmCa 

Figure 6 Vertical profile of temperature in the atmosphere. 

4>(k) <* k' 11 ^ 3 , (22) 

assuming isotropic and homogeneous turbulence, k represents spatial 
wavenumber in rad/m. Temperature is a "passive additive" (Tatarski, 1961) in 
the atmosphere, meaning that it does not affect the dynamics of the turbulence. 
Thus, the three-dimensional spectrum of temperature fluctuations also follows 
the Kolmogorov spectrum. Since Eq. (21) relates the index of refraction to 
temperature, the spectrum of refractive index fluctuations is also Kolmogorov 
(Tatarski, 1961) 

= 0.033 C‘(z) K' 11/3 , (23) 
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where C n 2 (z) is a constant of proportionality indicating the strength of the 
turbulence. In reality, this spectrum only applies to an intermediate range of 
wavenumbers known as the inertial subrange, shown in Fig. (7). The 
boundaries at the high and low spatial wavenumbers are the inner and outer 
scales, respectively, and are discussed below. 



Figure 7 Atmospheric spectrum of refractive index fluctuations showing inertial 
subrange and outer and inner scales. 


With this understanding of the effect of turbulence on a 
propagating E-field, Tatarski (1961) derived expressions for the statistical 
properties of the amplitude, phase, and irradiance of a spherical wave 
propagating through a turbulent medium as follows. Consider a scalar 
component of the vector wave equation Eq. (6) and assume a harmonic time 
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dependence for the E-field 


E(f,t) * u{?) e lul . ( 24 ) 

Perform the time derivatives and introduce the wavenumber k = 2 nJX and the 
index of refraction 

n = cfifi ( 25 ) 

(where c = velocity of light in vacuum, p = magnetic permeability of the 
medium, e = electric permittivity of the medium ) to get the scalar wave 
equation 

V 2 u + k 2 n z (t) u - 0. < 26 ) 

Expressing the index of refraction as 

n{f) - 1 + ^(7), (27) 

where |n,(7)| « 1, and using the Rytov method of smooth perturbations that 
employs u = exp(¥), ¥ = ¥ 0 + ¥, + .... and u = u 0 + u, + .... the wave 
equation becomes 

V 2 ?, + 2V¥ 0 -V¥ 1 + 2 k 2 n^(f) = 0. (28) 

This has a solution 

= dV> - < 29 > 
2nu 0 (r) J y \t - t '| 
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For a spherical wave 


u 0 (|?|) = Qex PWU , (O = constant) (30) 


and for small wavelengths X « <(,, where ^ (the inner scale) characterizes the 
size of the smallest fluctuations in the index of refraction, the Fresnel 
approximation in Cartesian coordinates reduces Eq. (29) to 




y'V*' 2 (x'+yh z 

2 zz' (z-z') _ 


dV'. 


(31) 


The statistics of the turbulence appear in 'P 1 and n, in terms of spectral 
expansions where n, contains the three-dimensional spatial frequency spectrum 
of refractive index fluctuations, <!>„(k). The Rytov approximation introduces the 
log amplitude fluctuation X 

X H In -A, (32) 

*0 


where A and A 0 are the amplitudes of the turbulence perturbed E-field u and 
the free space E-field u 0 


u = A e is 
u 0 ■ A 0 e ,s °. 

Then, after an extended series of manipulations (see Tatarski, 1961) that 
assume local homogeneity and isotropy, the variance of the log amplitude 


(33) 
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fluctuation X of a spherical wave propagating over a distance L becomes 


X 2 = An z k z k j^dz <D„(ic,z) sin 2 [— 


where k = lid'k, and k = spatial wavenumber (rad/m), which is used here 
instead of spatial frequency f (1/m) to be consistent with previous work. 
Assuming that the irradiance follows a log normal distribution (Tatarski, 1961), 
the normalized irradiance variance is a function of the log amplitude variance 

= exp( 4 X 2 ) - 1. ( 35 ) 


The assumption of isotropic, homogeneous turbulence gives a 
Kolmogorov spectrum of refractive index fluctuations of the form 
0„(k,z) = 0.033 C„ J (z) k- 11q (Tatarski, 1961). Assuming, for computational 
convenience, that this k' 11 ' 3 power law dependence holds over the entire range 
of spatial wavenumbers k, and that the turbulence strength C„ J (z) is uniform 
along the path, then integration of Eq. (34) gives the log amplitude variance for 
a spherical wave as 

3? = 0.124 C 2 n k 7 >* L"' e . (36) 

The normalized irradiance variance of Eq. (35) becomes 
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°L 

7 2 


(37) 


= exp(4 X 2 ) 1 = exp( 0.497 C 2 k 7/e Z. 11 ' 6 ) - 1. 


Low turbulence (small C n 2 ) and/or short path lengths make the exponent small 
enough to apply the approximation exp(a) « 1+a, giving 


7 2 


= 0.497 C 2 Ar 7/B L"* * p 2 e , 


(38) 


where the new parameter p D 2 , defined by this equation, serves as the baseline 
normalized irradiance variance for comparing modifications to the spectrum of 
refractive index fluctuations. p o 2 also facilitates plotting normalized irradiance 
variance over a broad range of integrated turbulence strengths. 

The above steps characterize the effects of the turbulent medium 
upon the irradiance statistics of a spherical wave assuming a simple 
Kolmogorov spectrum of refractive index fluctuations. Incorporating a high 
spatial frequency rolloff at the inner scale, as shown in Fig. (7), modifies the 
irradiance statistics. The inner scale ^ represents the physical size where 
viscosity of the medium smooths the velocity fluctuations by dissipating kinetic 
energy and thermal diffusion smooths the temperature fluctuations (thus 
refractive index fluctuations), removing the turbulent character of the medium at 
this scale. This smoothing causes a rolloff of the high spatial frequency energy 
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of the refractive index spectrum, Flatte, Wang, and Martin (1993) write the 
three-dimensional isotropic spectrum as 

= 0.033 C 2 (z) K' 11 ' 3 {39) 

where F(k{J represents a particular functional form of the inner scale. The 
parameter xt,, consists of the spatial wavenumber k (radians / meter) and the 
inner scale parameter (meters). For zero inner scale (i.e. no high spatial 
frequency rolloff), 

F(K{ 0 ) = 1, 0 < K < ~. (40) 

For theoretical and computational convenience, a Gaussian rolloff inner scale is 
often used (Tatarski, 1961) to represent the high spatial frequency rolloff from 
viscosity 

) ■ < 41) 

The more realistic viscous-convective enhancement inner scale advocated by 
Hill (1978) and Frehlich (1992) exhibits an enhanced spectrum for the 
wavenumbers slightly less than the rolloff point. Viscosity attenuates the kinetic 
energy and lowers the velocity fluctuations over regions near the inner scale 
size and smaller before thermal diffusion can smooth all the temperature 
fluctuations within these regions. This disparity alters the temperature and 
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corresponding index of refraction spectra. Inner-scale-sized patches of air have 
different temperatures but have little internal velocity variation. Since the 
Kolmogorov spectrum is based upon the spectrum of velocity fluctuations, these 
residual temperature fluctuations cause an enhancement to the Kolmogorov 
spectrum around the inner scale size. At higher spatial frequencies beyond the 
inner scale, thermal diffusion does smooth out these residual temperature 
fluctuations and the spectrum falls off sharply. 

Hill provides a plot of F(k^) versus for the viscous-convective 
enhancement, that is suitable for implementing the viscous-convective 
enhancement inner scale for numerical integrations and computer simulations 
(Flatte, Wang, and Martin, 1993). Frehlich presents a similar characterization of 
the viscous-convective enhancement inner scale based upon laser scintillation 
measurements and provides a four parameter fit to describe this version of the 
viscous-convective enhancement inner scale. 

Since most numerical simulations, such as the one used here, 
utilize discrete Fourier transforms, the spatial frequency grid mesh chosen for 
calculations introduces a maximum spatial frequency established by the 
Nyquist sampling criterion (discussed in Chapter III). This limit creates a grid 
cutoff inner scale at k^,,, 


F(k{ 0 ) = < 


(42) 
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Figure (8) illustrates these five inner scales for eight inner scale 
sizes. Solid curves represent the Hill version of the viscous-convective inner 
scale, dashed curves represent the Frehlich version of the viscous-convective 
inner scale, and dotted curves represent Gaussian inner scales. The solid line 
with the box shape represents the numerical grid cutoff inner scale at = 
318 rad/m for a 1024x1024 mesh. The inner scale values of 2, 3, 4, 5, 6, 7, 
10, and 15 cm refer to stratospheric propagation over 200 km, assuming an 
optical wavelength of 500 nm. 

The outer scale L 0 , corresponding to the spatial wavenumber 
K m = 2 tc/L 0 , represents the upper size limit to which the Kolmogorov spectrum 
applies. For scale sizes larger than L 0 (or spatial wavenumbers less than kJ, 
the refractive index fluctuations level off to a finite value as k approaches zero. 
Physically, such a limitation exists as the spatial wavenumbers approach zero 
because the refractive index fluctuations cannot become arbitrarily large, or 
equivalently, the energy represented by the spectrum must remain finite 
(Tatarski, 1961). Estimates of the outer scale for the stratosphere lie in the 
range of tens to hundreds of meters. However, the outer scale proves to be 
much more problematic to include in the spectrum of refractive index 
fluctuations because the atmosphere is anisotropic at these large sizes 
(Tatarski, 1961). The Von Karman spectrum (Tatarski, 1961), 


25 



F (Kl) 







(43) 


*„0O ■ 


0.033 C* 

(k 2 + 


-(f) 2 

e m . 


which has been used to incorporate an outer scale for some calculations, is 
based upon computational convenience more than physical understanding. 
Additionally the computer simulation results were compared against the Rytov- 
Tatarski predictions, which were derived without an outer scale. Therefore, 
these investigations did not incorporate an outer scale. 

2. Atmospheric MTF and Coherence Length 

Propagation through a turbulent medium not only affects the 
irradiance statistics, but also reduces the spatial coherence of the wave as 
characterized by the transverse coherence length. Fried (1966, p. 1372-1379) 
derived a long exposure modulation transfer function (MTF) for a spherical 
wave propagating through a turbulent atmosphere and used this to 
parameterize the E-field transverse coherence length, r 0 . Following the method 
and notation of Fried (which used spatial frequency f in ( m' 1 )), consider the 
spatial Fourier transform of the intensity of an E-field propagated through the 
turbulent medium and imaged by a thin lens 

t(/) = B f dx u’(x) u(x) e l2n!s , (44) 

where u(x) represents the E-field in the image plane and B represents a 
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normalization constant. Fried calls this the "MTF of the image-forming optical 
system", assuming a unit impulse (point) illumination and anticipating the fact 
that the ensemble average of t (f), which incorporates the effects of 
atmospheric turbulence, turns out to be real. Utilizing the Fourier transform 
property of a thin lens, 

u[x) = A f dv U{v) (45) 

where A is another normalization constant and U(v>) represents the E-field in 
the plane of the thin lens aperture, the MTF becomes 

t (?) = A 2 B f (tt U'{9 - XRf) U(v). ( 46 ) 

Expressing the E-field U(v) as the product of a zero turbulence propagation 
part, W(P) , and an atmosphere-induced perturbation, V(£) , 

U(v) = W{Q) V(v) = W{v) e'M * '*&. < 47 ) 

W( &) represents the uniformly illuminated aperture function. t( v) represents 
fluctuations in the log amplitude of the E-field and <(>(/) represents phase 
fluctuations. 

Taking the ensemble average < > of Eq. (46) over many 
realizations and substituting in Eq. (47) gives the long exposure MTF 
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<*(>)>*»* «**« = A * B / ^ w '(* - If#) W W < V*(* - Afl/) 1/(V) >.(48) 

Fried denotes the expectation of the fluctuations as the "atmosphere's MTF", 

MTF^JLf) = ( V'{v - \Rf) V(v) ), (49) 

and Hufnagel and Stanley (1964) call this the average mutual coherence factor. 
Substituting from Eq. (47), 

( V'(v - XRf) V(v) ) = * '(♦(»)-*(»-**f)l ), (50) 

Using the fact that 4> and t are Gaussian random variables, Fried shows that 
the atmosphere's MTF reduces to 

,51 > 


where D(XR |?|) represents the wave structure function and is related to the 
phase structure function and the log amplitude structure function D, by 

D(r) = D,(r) + D^(r), where r=\9-XRf\. (52) 

The log amplitude structure function D, and phase structure function are 
defined as 


D t (r ) - ( [|(P) - «(t >-XR?)] z ), 
D.(r) - < [*(0) - <M0 *fl?)] 2 >. 


(53) 
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and describe the expected variation in the log amplitude and phase, 
respectively, at two points a distance XRf apart. Denoting the MTF of the 
aperture function \N(v) as 

t 0 (/) ^ A 2 B [ dv I Y'(v-\R?) W{v), < 54 ) 


and using the fact that the atmospheric MTF of Eq. (51) is independent of v , 
the long exposure MTF of Eq. (48) now becomes 


^(h' ! lons exposure 


^f) e * 


DQ.R |f|) 


(55) 


or, 


A ™ exposure 


(56) 


The long exposure MTF is related to the mutual coherence of the 
E-field. The mutual coherence function (MCF) describes the autocorrelation 
between the E-field at two points and is defined as (Goodman, 1985) 

MCF ? ( U'Vvti) Wz’tJ >. (57) 

where U represents the complex E-field. For these investigations t,=t 2 and 
explicit time notation will be dropped. Assuming homogeneity, the MCF 
becomes the spatial autocorrelation of the E-field 
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MCF{f) = / eff' U'(F) U{f + ?). 


(58) 


Taking the ensemble average of Eq. (46) shows that the long exposure MTF 
comes from the average mutual coherence function (autocorrelation) of the 
E-field in the aperture of the lens 

<i(?)> = A Z B ( { dv U'(v - kR?) U(v) ) 3 < ( U'(v) U(v+kRt) )). (59) 
Rewriting Eq. (56), 

e \ = = < < U(F) U(?'+t) ) ) (60) 

< w(f) w(r'+r ))' 

Equation (60) is suitable for investigating the wave structure function D(r) since 
computer simulations provide the E-fields U and W on the right hand side of 
Eq. (60). 

Fried (1966, p.1380-1384) derived the wave structure function for 
a spherical wave based on the three-dimensional spatial frequency spectrum of 
refractive index fluctuations, <t n (K,z), 

D(r) = 8 n*k z dz / o “ [1 - jjZZ)] k (k, (61) 

where z represents the position along the optical path length, 0 < z < L . Fried 
analytically solves Eq. (61) for a spherical wave and simple Kolmogorov 
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turbulence that has an index of refraction structure function 


D„(r,z) = C 2 „(z) r 2 ' 3 , (62) 

and a three-dimensional spatial spectrum of refractive index fluctuations 

0(k,z) = 0.033 C„(z) iT 11/3 . (**) 

Substituting Eq. (63) into Eq. (61) and integrating over spatial frequency k gives 
D(r) = 2.91 k 2 r 5 ' 3 [ q L dz Cfc) (|f /3 . (64) 


Assuming C„ J (z) = constant along the optical path gives 

D(r) = 1.089 k 2 C 2 n L r sl3 . ( 65 ) 


This equation relates the wave structure function D(r) to physical parameters 
k=2n/A., C n J , and L for the propagation through turbulence. Following Fried and 
defining the constant r c 


6.88 | 3 ' 5 

1.089 k 2 C 2 „ U 


( 66 ) 


Eq. (65) becomes 

D(r) = 6.88 (-^) 5/3 . (67) 

r o 

The parameter r Q characterizes the coherence length of the E-field because the 
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resolution allowed by the turbulent atmosphere is Q R * A. / r 0 (Fried, 1966, 
p. 1380-1384). 

Substituting Eq. (67) into Eq. (60) gives 

0 ~ 3 44( ? )8/3 = ( ( U(t>) UV'+I) ) ) (68) 

( Wit') W{t'+f) > ’ 

and now directly relates a measure of the coherence length, r 0 , to the E-fields 
produced by computer simulation. Similarly, substituting the integral form of the 
wave structure function, Eq. (61), into Eq. (60) and numerically integrating for 
spectra <l> n (ic,z) with different inner scales allows comparing theory and 
computer simulation. This comparison facilitates validation of the computer 
simulations incorporating an inner scale at low turbulence strengths where the 
perturbation-based theory remains valid, and provides a mechanism to explore 
conditions of high turbulence strength and/or long propagation path lengths (the 
saturation regime) where the perturbation theory may no longer hold. 
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III. COMPUTER SIMULATION 


A. PROPAGATION CODE 

These investigations utilized a wave optics computer simulation code 
written by Brent Ellerbroek at the United States Air Force Phillips Laboratory in 
Albuquerque, New Mexico that incorporates phase screen generation routines 
written by Greg Cochrane, also with the Phillips Laboratory. The code, known 
as YAPS (Yet Another Propagation Simulation), is a general purpose adaptive 
optics simulation code written in FORTRAN that models optical propagation 
through a turbulent atmosphere, sensing of the wavefront with Hartmann 
elements (Hudgin, 1977) and a CCD array, optimized phase correction 
calculation, and phase compensation via a deformable mirror, all within a time 
indexed framework. These investigations utilized only the propagation portions 
of the code, modified the code to parameterize turbulence strength with C„ J 
instead of r 0 , added different source configurations, incorporated Gaussian, Hill, 
and Frehlich inner scales, and reduced the memory requirements. 

The computer simulation utilizes the split-step method to simulate 
propagation of the E-field through turbulence, as illustrated in Fig. (9). A source 
E-field (left) was propagated in steps (with zero turbulence) over the 
propagation distance L with a random phase screen applied to the field at each 
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Figure 9 Computer simulation using the split-step method. The source (left) is 
propagated in steps out to z = L, with phase screens applied at each step. 


step. This method was based upon the extended Huygens-Fresnel principle 
(Yura, 1992) 

E p = y' /T E a K(Q) o'* dA, (69) 

which is Eq. (10) with an extra e* factor that incorporates the random variations 
in log amplitude { and phase <|> 

s e 1 e 1 *. (70) 

For paraxial propagation, if the distances n are short enough, diffraction and 
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interference will not have a chance to cause significant log amplitude 
fluctuations t, allowing the approximation 

gf * g> ♦. (71) 

Inserting this approximation into Eqs. (11) - (19) gives 

E{t,z) = fdfee~'* kz Ifl* /dp E(p,0) *'♦»> ©-' 2 ”'p. (72) 

Thus, a single step of the split-step method requires applying the phase screen 
e 1 * to the E-field and propagating the result a distance z using the Huygens- 
Fresnel propagator. 

To implement this sequence of steps, the YAPS code followed a user- 
defined list of tasks (stored as an input file) and called the appropriate routines 
to accomplish those tasks. For these investigations, a typical list (see the 
Appendix) first set up the initial parameters, including random number seed, 
wavelength, number and size of fields, source characteristics, and phase 
screen size. The list then initialized the field grid and applied the initial source 
distribution, propagated the field in steps, generated and applied phase screens 
between steps, and finally saved the complex values of the propagated field. 

The YAPS propagations were run on Sun SparcStation-10's having 128 
megabytes of RAM. This memory constraint allowed a maximum grid size of 
1024x1024 (choosing N as an integer power of 2) for the current version of 
YAPS. Each run consisted of 30 propagations, each with 32 phase 
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screens/steps, and required -16 hours to complete. These investigations 
normally utilized nine SparcStation-10's simultaneously doing runs with different 
parameters. In all, these investigations consumed roughly 6000 hours of 
computer time. A small Cray mainframe was available with more RAM, but was 
not utilized extensively because the multiple SparcStations provided more 
overall computational capacity. 

Once the YAPS code had generated and saved the realizations of the 
E-field propagated through turbulence, separate routines written in Interactive 
Data Language (IDL) analyzed and displayed the fields. Analysis included 
displaying two- and three-dimensional plots of the intensity field, calculating the 
dependence of intensity and normalized irradiance variance on radial distance, 
calculating the normalized irradiance variance over the central portion of the 
field, calculating the atmospheric MTF and corresponding coherence length, 
and calculating Strehl and intensity ratios (defined in Section E below). 

The wave optics computer simulation required many choices to model 
the stratospheric propagation scenario and to ensure validity of the simulation. 
These choices included wavelength, propagation distance, inner scale, grid size 
N, the physical size of each grid element Ax, the maximum strength of 
turbulence p 0 2 max , the source function of the E-field, the number of phase 
screens, low spatial frequency corrections, number of realizations, and methods 
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of calculating average statistical values. The following sections address these 
and other choices and develop guidelines for validity of the simulations. 

B. PHYSICAL PARAMETERS 

Wave optics computer simulation of the E-field requires selection of the 
parameters that describe the physical properties of the propagation. Usually a 
specific propagation scenario has been selected for modeling, giving the 
desired wavelength X, wavenumber k=2nlX, propagation distance L, and the 
range of turbulence strengths C„ 2 . The stratospheric propagation scenario 
chosen here used X = 500 nm, L = 200 km, and C„ 2 = [1x1 O' 21 , IxlO' 16 ] m" 2 ' 3 . 
From these physical parameters, other useful scaling parameters occur, such 
as the Fresnel wavenumber (Martin and Flatte, 1988) 

" ■ 1 If)"’ ■ jin:’ |73) 

the Rytov-Tatarski normalized irradiance variance (Flatte, Wang, and Martin, 
1993) 

^ - a Cl k 7 " = p* , (74) 

where a=1.23 for plane waves and a=0.497 for spherical waves, and Fried’s 
coherence length for spherical wave propagation (Fried, 1966, p.1380-1384) 
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r ,(_5*5_f S 

(?5) 

= (- 6-38 ) 3/5 (for constant C 2 ). 

1.09 k 2 C„ U 

The inner scale of the turbulence, «<,, often is not known exactly, but values can 
be estimated. For stratospheric propagation, estimates of the inner scale are 
around 1 - 15 cm (Beland, 1993). Once the inner scale is known, the Hill and 
Frehlich versions of the viscous-convective enhancement inner scales (Hill and 
Clifford, 1978) (Frehlich, 1992) were implemented with the parameter , 
where k represents spatial frequency. The outer scale L 0 for stratospheric 
propagation lies in the range of tens to hundreds of meters but, again, due to 
the difficulty in parameterization, was not included in these simulations. 

C. COMPUTATION GRID 

The E-field was represented as an NxN array of complex numbers in the 
plane perpendicular to the axis of propagation . Generally, N was chosen as 
large as possible, consistent with the amount of computer memory available 
and the time required to run a simulation. Larger grid size N allowed a wider 
spatial extent of the field and/or sampling of higher spatial frequencies (denser 
mesh). The physical size that each grid element represented had to be chosen 
to ensure validity of the computer simulation. 
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Knepp (1983) discussed the relation of grid element size Ax and physical 
grid extent NAx to inner and outer scales, to proper sampling of the phase 
screens, to angular spreading of the field, and to proper sampling of the spatial 
frequency quadratic phase factor in the Fourier transform formulation of the 
Huygens-Fresnel propagation. For the latter, the quadratic phase factor in 
Eq. (20) takes the form 


4 > 


k 2 L 
2 k’ 


(76) 


where k is a spatial wavenumber in rad/m. Applying the Nyquist criterion, 
which requires that this quadratic phase change by less than 71 across one grid 
element Ax, Knepp derives 

L < 2N[Axf (77) 

A. 


Roberts (1986) applied similar techniques and derived Eq. (77) without 
the factor of 2. Again applying the Nyquist criterion to the quadratic phase 
factor in the Fourier transform propagator Eq. (20) but utilizing more succinct 
differential notation, 
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and, 


A(J> * Af2 rt A L < n, 


(79) 


where f max represents the spatial frequency with the maximum rate of change of 
phase. At the edge of the grid, f^, is 


1 

2 Ax' 


(80) 


After substituting and rearranging, 

L < N (Ax) 2 

A, 


(81) 


When the propagation distance L is Known, as for a specific propagation 
scenario, and when the grid dimension, N, is known, then the grid element size 
is 


Ax > 


no. 

i n 


(82) 


Sampling of the phase factor in the spatial frequency domain must meet 
the Nyquist criterion as a minimum. Some circumstances may warrant applying 
even stricter sampling criteria such as restricting the phase to change by less 
than an across one grid element, 0 < a < 1. Furthermore, suppose the 
propagation involves spatial frequencies out to pf ma „ where 0 < p < 1. Then 
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A4> < an, and the maximum spatial frequency is pf max . Substituting these 
constraints into the analysis gives 


Ax 


[HI 

i N a 


(83) 


Whichever is used, Eq. (82) and (83) provide simple formulas for choosing the 
minimum grid element size for the propagation. 

Spherical wave propagation places an additional constraint on the choice 
of grid element size. In the parabolic approximation, a spherical wave has a 
quadratic phase curvature which represents divergence from a point source a 
distance S (focal distance) away 


where p measures the radial distance from the propagation axis. Roberts (1986) 
analyzed the sampling criteria for this phase just as for the spatial frequency 
phase (though applied to the case of a two-step Fourier transform Huygens- 
Fresnel propagator and not applied specifically to spherical waves). Proceeding 
as above and applying the Nyquist criterion to this phase factor 

_ _d_ - Ai, (85) 

dp dp i S Ap 
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< n. 


(86) 


A(j> 


A P 2 * Pmax 
X S 


where p max represents the maximum radial distance. Using p maK corresponding 
to the edge of the grid and Ap = Ax, 


N Ax 
2 


(87) 


Substituting Eq. (87) into Eq. (86) and rearranging, 


Ax < 




(88) 


Again generalize the analysis by requiring the maximum phase change 
across one element to be an (0 < a < 1), and requiring the field energy to be 
confined within a region of radius yp max , (0 < y < 1). This gives 


Ax 


[HI 

N N y- 


(89) 


Combining Eq. (83) and (89), 


[HI 

N N a 


< Ax < 


I X Sa 

N N y' 


(90) 


Spatial frequency sampling considerations for the Huygens-Fresnel propagator 
have placed a lower limit on Ax, while spatial sampling considerations of the 
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quadratic approximation for spherical wave phase have placed an upper limit on 
Ax. 


The computer simulations used in these investigations examined high 
levels of turbulence that introduced energy into the highest spatial frequencies 
representable on the grid and that scattered energy to the edges of the spatial 
grid. The simulations also assumed that the Nyquist sampling criterion was 
sufficient. These considerations specified the parameters a = p = y = 1, and 
lead to 




Additionally, these simulations propagated a spherical wave from the point 
source at the origin (focus) out to a distance S, i.e. S = L, making the 
inequalities in Eq. (91) become the equality 




These choices unambiguously determine the guideline for grid element size for 
spherical wave propagation based upon sampling considerations in the spatial 
and spatial frequency domains. Equation (92) also determines the grid element 
size for plane wave propagation since minimizing the grid element size of 


44 



Eq. (82) optimizes the sampling of high spatial frequency distortions caused by 
turbulence. 

Some simplifications were made in the above analysis. First, the 
maximum spatial distance and maximum spatial frequency used in the 
quadratic phase factors were chosen for the nearest edge of the grid and 
correspond to the radius of the largest circle that will fit inside the grid. These 
choices disregarded the corners of the grid, but this omission should not have 
affected the simulations greatly because the majority of the energy in both the 
spatial and spatial frequency domains was confined within the radius of the 
circle to minimize aliasing, Second, the Nyquist criterion was applied to the 
phase change across the x or y dimension of the element. Analyzing the phase 
change between opposite corners of a grid element increases Af or Ap by . 


Finally, an NxN grid with N even requires placing the (x,y) = (0,0) point at a grid 
point, such as (N/2 +1, N/2 +1) that is not the exact geometric center of the 
grid. The distance to the nearest side is (N-2)/2 elements instead of N/2, which 
is a minor change for large N. Incorporating these three considerations into the 
analysis for a spherical wave gives 


2 \ L (A/-2 ) < < 

\ N 2 


N 2 (A/-2) 


(93) 
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However, squaring Eq. (93) and rearranging leads to 

L < 1 - S, (94) 

4 (A/-2) z 

so that a spherical wave propagation with these additional constraints should 
only be carried out over a maximum distance less than S/4. Since these 
investigations needed to propagate a spherical wave over the full focus 
distance S starting at the source, the simpler expression Eq. (92) was used to 
determine the grid element size for these investigations. As a comparison, 
Flatte, Wang, and Martin (1993) used a grid element size of 0.7mm for a 
1024x1024 grid and 0.5mm for a 2048x2048 grid with XL = 0.000638 m 2 . 
Equation (92) with these parameters prescribes grid element sizes of 0.78 mm 
and 0.55 mm, respectively, which are ~10% larger to meet sampling 
considerations for the Huygens-Fresnel propagator. 

The split-step method in the computer simulation divides the optical path 
into multiple steps. However, the distance L used to determine the grid 
element size must correspond to the total propagation path length and not the 
step size. To justify this, consider the vacuum propagation of a field across the 
distance L = Az. If propagation occurs in a single step, then the Fourier 
transform Huygens-Fresnel propagator, Eq. (20), becomes 

E(Az) - !FT[ e ET[E(0)]]. < 95 > 
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If the path has two steps, then 


E(Az) = IFT[ e 


FT[E(Az/2)] ], 


and 


E(Az/ 2) = IFT[ e 2 F7[E(0)]]. 


( 97 ) 


Substituting Eq, (97) into Eq. (96), 


/nX 

E(Az) = IFT[ e 2 FT[ IFT[ e 2 FT[E( 0)]]]]. 


(98) 


But the middle Fourier transform/ inverse Fourier transform pair cancel each 
other, and this two step vacuum propagation becomes the one step propagation 
of Eq. (95). Thus, the grid element size must be chosen to correspond to the 
total distance Az = L since multi-step vacuum propagations mathematically 
reduce to a single step of L. 


D. SOURCES 

The choice of a given propagation scenario significantly affects the 
irradiance and coherence statistics. Plane wave propagation differs from 
spherical wave propagation, as seen in the Rytov-Tatarski theory (Tatarski, 
1961), while beam wave propagation (Gaussian intensity profile) exhibits an 
intermediate behavior (Ishimaru, 1978). The source must also be chosen to 
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keep the lateral extent of the source propagation within the spatial and spatial 
frequency limits of the computer simulation. 

Coherent plane, beam, and spherical waves differ in the shape of the 
isophase surfaces of the E-field. Beam wave sources possess a quadratic 
phase (Ishimaru, 1978), 

<Mp) = P Z . 09) 

where R 0 represents the radius of curvature of the wavefront. Plane wave 
sources have an infinite radius of curvature so that 

<|>(p) = constant. ( 10 °) 

Spherical wave sources have a radius of curvature equal to the propagation 
distance from the origin, or focus, S 

■mp) - mod 

Turbulence introduces random phase shifts across the wavefront while 
diffraction and interference further distort the wavefront during propagation 
causing the light to become partially coherent. The resulting partially coherent 
E-field no longer possesses a simple isophase surface and the coherent 
wavefront characterizations no longer apply. The differences between plane, 
beam, and spherical propagation appear in differences in statistical properties, 
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such as normalized irradiance variance (Tatarski, 1961) 


2 1.23 Cl k 7/t L 11 *, plane 

^ = { ( 102 > 

7 0.497 Cl k 71 * L 11/6 , spherical 

(see Ishimaru, 1978 for the more complicated beam wave expression), and the 
coherence length (Fried, 1966, p. 1380-1384) 


3.02 [Cl) W k^ 5 L 3 * 5 , spherical 
' 1.68 {Clr 315 kL-W, plane. 


When the E-field has propagated through turbulence to the far field of 
the source, diffraction has caused the E-field to expand laterally so that its 
statistical properties approach those of a spherical wave. The far field begins 
at a distance (Saleh and Teich, 1991) 

z « 10 (104) 


where p is the radius of the source. Spherical propagation properties result 
when the maximum radius of the source approximately satisfies 

Pmax < y'O.I X z. (1° 5 ) 

For a stratospheric propagation scenario with L=200 km and X = 500 nm, the 


49 



source radius must be less than 10 cm to obtain enough divergence to produce 
spherical propagation statistics at z = L. 

Plane wave propagation statistical properties result when the E-field is 
evaluated in the near field of the source because the E-field does not have the 
opportunity to diffract or expand significantly. Starting with Eq. (105) for the far 
field point, near field propagation satisfies 

Z „ J_ « 10-H! (106) 

10 X X 


or, 

p™ > v'WTT < 107 > 

For the stratospheric propagation scenario, source radius p m(K > 1 m for near 
field propagation. However, validity of the Fresnel approximation to the 
Huygens-Fresnel equation places an upper bound on (Saleh and Teich, 
1991) 

(pL) 2 « 4 pr 3 X, (108) 

(derived from considering the Taylor series expansion of the Cartesian 
expression for r). For the L = 200 km, X = 500 nm stratospheric propagation 
scenario, the upper bound on source lateral extent becomes p^, < 350 m. A 
1024x1024 grid with the grid element size Ax = 0.99 cm from Eq (92) gives a 
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maximum grid radius of only 5 m, so that the Fresnel approximation is satisfied 
for any source represented on this grid. 

Ishimaru (1978) evaluated the irradiance statistics for the intermediate 
Gaussian beam wave case. The normalized irradiance variance along the 
beam axis depends upon the beam waist size, W 0 , and the radius of curvature, 
R 0 . Numerical integration of the log amplitude variance formula (Walters, 1994) 
provides a smooth transition from spherical wave variance statistics to plane 
wave statistics, with a dip in between, as shown in Fig, (10) . Numerical 
simulation results with Gaussian and Airy-type sources of varying widths have a 
corresponding behavior, as shown in Fig. (11). 

The finite grid in a computer simulation places limitations on the source 
E-field. The physical grid element size sets a lower bound on the width of a 
narrow source approximating a point source. A source of width close to a 
single grid element may still be undersampled, causing the resulting propagated 
field to exhibit sidelobes from absence of the proper high spatial frequencies in 
the representation. Correspondingly, a spherical wave from a very narrow 
source will propagate with large angular divergence that may exceed the 
physical dimensions of the grid and cause energy to leak off the grid and be 
aliased back into the field. Thus, the source must be wide enough to constrain 
the propagated field so that most (>90%) of the energy remains within the grid. 
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Figure 11 Normalized irradiance variance versus beam waist size W 0 , where 
Q=tcW 0 2 /A.L. Solid line = Airy-type source; dashed line = Gaussian source. 





Similarly, a plane wave source cannot extend too near the grid edge because 
turbulence can scatter energy off the grid only to be aliased back in. 

Statistical calculations on the final beam pattern require a reasonably 
uniform central patch to perform computations. If the region over which 
statistical calculations are made has variations, such as sidelobes arising simply 
from vacuum propagation, these variations become a part of the statistics, such 
as the normalized irradiance variance. Often, Gaussian type profiles present a 
minimum of such variations because they do not have as much energy in high 
spatial frequencies. However, because they are not flat over the calculation 
region, different regions of the beam must often be weighted with respect to 
their mean intensity in doing statistical calculations such as the normalized 
irradiance variance. 

To meet these constraints, Martin and Flatt6 (1990) applied a quadratic 
phase curvature to a Gaussian source, thereby increasing the divergence and 
flattening the final irradiance pattern. Spherical wave sources were simply 
chosen narrower than the plane wave sources to make them diverge more. 
Flatte, Wang, and Martin (1993) also used a super-Gaussian to model an 
extended beam source. 

The simulations presented here used an alternate method suggested by 
Ellerbroek (1993). The final field at z = L was specified as an aperture of 
radius equal to one half the grid radius and given a quadratic phase 
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corresponding to a spherical wave originating at the origin z = 0 (Fig. (12)). 

This aperture field was then backward propagated without turbulence from z = L 
to the focus z = 0, effectively Fourier transforming the field and yielding the 
numerical equivalent of the Airy pattern (Fig. (13)). This source was then used 
as the approximation to a point source for all spherical wave propagations. The 
advantages of this method were that the source was quite narrow, having 
appreciable amplitude over only a few grid elements, and that the zero 
turbulence propagated field at z = L was necessarily constrained within the grid 
and possessed a central region with uniform illumination. 

The significant energy at high spatial frequencies required to represent 
the sharp edges of the initial aperture at z = L created one drawback since it 
caused strong Fresnel fluctuations at intermediate distances between 0 and L. 
To investigate the significance of these Fresnel fluctuations, the energy at these 
high spatial frequencies were reduced by windowing the aperture before the 
backward propagation by applying a Gaussian rolloff to the cylinder edges, 

( 1 , r< 0.7 r' 


1*1 = 


-]") , r*0.7r', n= 2, 


where r' represented the corresponding aperture radius. Because this rolloff 
reduced the energy in high spatial frequencies for the source representation, 
this slightly increased the size of the final central region of the E-field where 
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Figure 12 Cross section of intensity profile for final propagated beam showing 
confinement of energy and uniform central illumination (beam radius = 70 grid 
elements here). 



Figure 13 Cross section of intensity profile of numerical equivalent to Airy 
source. 
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statistical calculations could be performed. Other values for the power n on the 
exponent were tried (n = 1, n = 5, n = 10) but the simple Gaussian (n = 2) 
worked best. While this rolloff appeared useful, identical runs, one with the 
sharply defined aperture edge and one with a rolled-off Gaussian edge, showed 
less than 0.5% difference in the normalized irradiance variance values. This 
difference was negligible compared to statistical fluctuations in normalized 
irradiance variance of up to ~10% between runs with different random phase 
screens. Empirical observation of simulations revealed that once even modest 
amounts of turbulence existed (typically (3 0 2 > 0.05 and well within the Rytov 
regime), the high spatial frequency energy introduced by the turbulence 
dominated the details of the source. 

Another deficiency was the representation of the initial cylindrical 
aperture on a rectangular grid. Due to the Cartesian nature of the grid, the 
curved beam edge was actually jagged instead of circular. However, for large 
enough grid (e.g. 1024x1024) this departure from a cylinder introduced 
negligible effect. The only case where a difference was noted was with a 
256x256 grid, and increasing the beam radius a small fraction of a grid element 
size eliminated that difference. 

Another limitation of this method was that the final E-field could not 
approach the grid radius. If this happened, then high turbulence strengths 
would scatter energy off the grid, producing aliasing. Due to the periodic 
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Fourier transform method of propagation, this energy would not be lost, but 
aliased back in at lower spatial frequencies. Martin and Flatt§ (1990) 
implemented an attenuating region just inside the grid radius to absorb this 
energy and to prevent its being aliased back in. This obviously introduced a 
type of windowing function in the spatial domain and did not conserve energy in 
the field, possibly complicating final field comparisons. Rather than include 
those effects, these simulations simply chose the aperture radius sufficiently 
small (at one half the grid radius) to prevent significant energy scattering off the 
grid in the spatial domain while still providing a relatively uniformly illuminated 
central region for statistical calculations. 

Obviously, other choices for a source existed. For example, the initial 
field could have been expressed according to an analytic expression for the 
Airy pattern (Fig. (14)). This source gave a vacuum propagated final field that 
was very close to a broad cylinder, but which exhibited noticeable sidelobes at 
the edges (Fig. (15)). These sidelobes came from spatial truncation of the Airy 
pattern. The back propagated numerical Airy pattern of Fig. (13) differs from an 
analytic Airy pattern in a way that eliminates the sidelobes in the final irradiance 
field. 

A point source could be modeled by a single, nonzero pixel at the center 
of the grid being nonzero, as shown in Fig. (16). Figure (17) shows that the 
vacuum propagated field exhibited noticeable sidelobes at the edges and 
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Figure 14 Cross section of intensity profile of analytic Airy source. 



Figure 15 Cross section of intensity profile from analytic Airy source after 
vacuum propagation. 
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Figure 16 Cross section of intensity profile for the single-grid-element point 
source. 



Figure 17 Cross section of intensity profile for single-grid-element point source 
after vacuum propagation. 
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covered the entire aperture, meaning that turbulence would immediately scatter 
energy out of the grid only to be aliased back in and distort the field. A similar 
source was a uniformly illuminated grid at z = L with quadratic phase, backward 
propagated to z = 0. While this eliminated the vacuum propagation ringing 
problem by definition, it immediately suffered from aliasing when turbulence was 
nonzero and scattered energy off the grid. 

E. MAXIMUM TURBULENCE STRENGTH 

One of the more difficult choices in a computer simulation of propagation 
is determining the range of turbulence strengths over which a simulation is 
valid. In general, the processes of optical propagation through a turbulent 
medium are not band limited in spatial frequency. But when a computer 
simulation uses a finite grid to implement a source function, to propagate via 
the Huygens-Fresnel principle, and to model the turbulence by phase screens, 
aliasing inevitably occurs due to the finite sampling interval. This problem only 
becomes more severe as turbulence strength increases. The coherence length 
measures the physical distance over which the mutual coherence of the E-field, 

{E* E ), declines to e' 1 of its peak value. As turbulence increases, the E-field 
fluctuates significantly over smaller and smaller distances and the coherence 
length decreases. An NxN discrete grid representation of the E-field samples at 
a specific minimum distance, and hence at a maximum spatial frequency. If the 
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E-field fluctuates in less than this minimum distance, the grid samples will not 
represent the E-field accurately, causing aliasing. The question is not whether 
aliasing occurs, but how much aliasing occurs for a given set of simulation 
parameters, and at what point does aliasing invalidate the results of the 
simulation. 

These investigations identified three telltale signs of aliasing. First, as 
aliasing became significant, irradiance plots of the turbulent E-field lost structure 
and exhibited an isotropic fine-grained appearance, as Fig. (18) shows. 

Second, the irradiance that should have been roughly radially symmetric 


Figure 18 Intensity plot showing fine-grained pattern due to significant aliasing. 
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became bounded within a fuzzy rectangular region, as Fig. (19) illustrates. The 
irradiance at the center of the image increased, or started peaking, while the 
irradiance at the edges decreased, as shown in Fig. (20). 

Figure (21) also illustrates this peaking behavior by plotting average 
irradiance versus radius, in units of grid element size. The E-field propagated 
with zero turbulence had the radially symmetric aperture profile in this 
simulation. As the turbulence became stronger, the sharp edges of the initial 
cylindrical irradiance pattern became rounded and energy spread outward. 
However, when significant aliasing started occurring around (3 0 2 - 5 for this 
64x64 grid representation, the energy began creeping inward and the center of 
the field started increasing in irradiance. 

The irradiance at the edge of the grid was -1/40 that of the center when 
significant aliasing began. This indicates that very little energy leaked off the 
spatial grid due to beam divergence. Actually, the use of Fourier transforms in 
the propagation preserved the energy so that any energy that leaked off the 
grid was aliased back into the field on the grid. Energy that leaks off in the 
spatial domain results from undersampling in the spatial frequency domain, and 
vice versa. 

Aliasing most often occurs from undersampling in the spatial domain, 
which corresponds to energy leaking off the grid in the spatial frequency 
domain. Prudent choice of source and the final irradiance patterns across the 
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Figure 19 Intensity plot showing boxed perimeter due to significant aliasing. 
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Figure 20 Intensity plot showing peaking of energy toward center due to 
significant aliasing. 
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Radial distance (bin f) 


Figure 21 Radial intensity profiles showing peaking behavior: p„ 2 = 0 (solid 
line), 0.5 (pluses), 1.5 (asterisks), 5 (dots), 15 (diamonds), 50 (triangles), 
150 (squares), 500 (X’s). 
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grid can minimize the amount of energy that leaks off the grid to be aliased 
back in for the spatial domain. However, as turbulence strength increases, the 
coherence length decreases and more energy occurs at higher and higher 
spatial frequencies. Eventually spatial undersampling occurs and this high 
spatial frequency energy leaks off the spatial frequency grid and is aliased back 
in. Figure (22) gives the radial power spectral density corresponding to Fig. 

(21) and shows the spread of energy to higher spatial wavenumbers as 
turbulence strength increases. The low turbulence spectra possess a strong 
central peak that becomes more rounded and flatter as turbulence strength 
increases, causing more energy to leak off the grid. Eventually, enough energy 
has leaked off the grid and aliased back in to make the spectrum approximately 
uniform at all spatial frequencies. 

Figure (23) actually shows this energy being reflected back in after it 
spills off. To reveal this phenomenon, the radial power spectrum for a 512x512 
grid has been divided by the corresponding part of the power spectrum on a 
1024x1024 grid. Since the larger grid has been chosen with a finer mesh than 
the 512x512 grid, the larger grid will not experience significant aliasing as soon. 
As turbulence strength increases, the ratio shows the enhancement of energy 
at the edge of the 512x512 grid (bins 200-256) as the energy leaks off / reflects 
back in. At very high levels of turbulence, the aliased energy has spread 
inward over all spatial frequencies. 
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Figure 22 Power spectrum radial profile as turbulence strength increases: 3 
= 0 (solid line), 0.5 (pluses), 1.5 (asterisks), 5 (dots), 15 (diamonds), 50 
(triangles), 150 (squares), 500 (X's). 










Aliasing appears to cause the intensity peaking behavior of the central 
field. These investigations parameterized this onset of peaking with five 
methods, which provided guidelines for maximum turbulence strength p 0 2 max 
valid for a given grid size. An irradiance ratio and a Strehl ratio were first used 
to identify the onset of the peaking. A Gaussian source was propagated at 
turbulence strengths P 0 2 = [5x10‘ 4 , 5x1 O' 3 , 5x1 O' 2 , 0.15, 0.5, 1.5, 5, 15, 50, 
150, 500] using 64x64, 128x128, 256x256, 512x512, and 1024x1024 grids. 

Ten runs for each case were used to calculate the average irradiance as a 
function of radius. These irradiance profiles were used to calculate an 
irradiance ratio, defined as the irradiance at the grid center divided by the 
irradiance at the maximum grid radius, and also the Strehl ratio, defined as the 
irradiance at the grid center for the E-field propagated through turbulence 
divided by the center irradiance for an E-field propagated through zero 
turbulence. Figure (24) plots the irradiance ratio versus turbulence strength for 
different grid sizes, and Fig. (25) shows the Strehl ratio versus turbulence 
strength for different grid sizes. 

As turbulence strength increased, diffraction and scattering spread the 
energy outward and caused the irradiance and Strehl ratios to decrease. 
Eventually, both ratios reached a minimum and then started increasing because 
of the alias-induced intensity peaking. This minimum occurred at higher 
turbulence strengths for larger grid size because the larger grid sizes sampled 
with a finer mesh and did not significantly alias as soon. Using estimates of 
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these minima as a guide to the onset of significant aliasing, Fig, (26) plots 
these p 0 2 values versus grid size from the irradiance and Strehl ratios, along 
with least squares fits and extrapolations to larger grid sizes. 

Martin and Flatte (1988) and others have validated computer simulations 
for predicting statistical properties such as normalized irradiance variance from 
the E-field propagated through turbulence. The departure of the computer 
simulation E-fields and intensities from their known or expected smooth 
behavior as turbulence increases represents another telltale sign of significant 
aliasing. Figure (27) shows the normalized irradiance variance calculated from 
computer simulated E-fields versus p o 2 using grids of sizes 64x64, 128x128, 
256x256, 512x512, and 1024x1024. The irradiance value used for 
normalization was taken as the average irradiance over the central calculation 
region. Though this single value normalization is not optimal (see the following 
section, F. Additional Simulation Parameters), it does make the normalized 
irradiance variance calculation sensitive to the peaking behavior because 
peaking adds large amounts of unphysical variance. In the Rytov regime 
(p 0 2 < 1.0), all grids successfully simulated the E-field as indicated by their 
normalized irradiance variances agreeing with the Rytov-Tatarski theory. In the 
saturation regime (p 0 2 > 1.0) the normalized irradiance variance saturates and 
turns downward, eventually approaching unity according to asymptotic theory. 
Assuming that the 1024x1024 grid provides the most accurate propagation 
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Figure 26 Irradiance ratio (diamonds / pluses) and Strehl ratio (triangles / 
asterisks) estimates and least squares fits showing maximum turbulence 
strength versus grid size for valid simulations. 




Beta squared 


Figure 27 Normalized irradiance variance versus turbulence strength from 
simulations with different grid sizes. 


simulation, the 64x64 grid normalized irradiance variance departs from the 
1024x1024 values before the peak was reached, and 128x128 grid normalized 
irradiance variance departed just beyond the peak. Their meshes were not 
small enough to adequately sample the E-field, and the resulting aliasing 
caused peaking that drove up the normalized irradiance variance calculation. 
The 256x256 grid produced the saturation peak, but soon suffered from 
significant aliasing. The 512x512 calculations successfully produced the peak 
but departed from the 1024x1024 predictions before p 0 2 - 50 . Judging from 
these behaviors, the 1024x1024 grid probably remains valid through p 0 2 ~ 50 . 
All grids eventually showed the normalized irradiance variance anomalously 
rising for high enough turbulence strength because aliasing produced peaking. 
Estimates of the maximum p 0 2 values at which the five grid sizes remained valid 
are plotted, along with a least squares fit extrapolated to larger grid sizes, in 
Fig. (28). 

In a similar manner, the coherence length and the half width at half 
maximum (HWHM) of the atmospheric MTF predicted using simulations with 
different grid sizes also show unphysical behavior as aliasing became 
significant. Figure (29) plots the simulation derived coherence length, 
normalized by the Kolmogorov turbulence theoretical values, versus turbulence 
strength measured by p 0 2 . Figure (30) shows the corresponding plot for 
HWHM. Both show the eventual strong rise in coherence length and HWHM 
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Figure 28 Estimates from normalized irradiance variance of maximum 
turbulence strength versus grid size for valid simulations. 
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Figure 29 Spherical wave coherence length from simulations versus turbulence 
strength for grid sizes 128x128 (diamonds), 256x256 (triangles), 512x512 
(squares), and 1024x1024 (X's). 
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Figure 30 Spherical wave HWHM from simulations versus turbulence strength 
for grid sizes 128x128 (diamonds), 256x256 (triangles), 512x512 (squares), and 
1024x1024 (X's). 




compared to theory for strong turbulence when aliasing produces peaking. 
Estimates of the minima of the plots can be used as another indicator of the 
turbulence strength at which significant aliasing occurs. Figure (31) plots these 
points versus grid size, along with a least squares estimate extended to larger 
grid sizes. 

The onset of peaking and departure from expected physical behavior 
provide symptoms of aliasing but do not indicate the amount of aliasing 
required to cause them. Determination of the fraction of total field energy that 
is aliased serves as one parameterization of the amount of aliasing and 
accomplishes two goals: (1) relate the amount of aliasing occurring to an easily 
measurable characteristic of the simulation , and (2) determine how much 
aliasing must occur to invalidate the computer simulation. 

First, to parameterize the amount of energy aliased, the same size 
Gaussian source was applied to 64x64, 128x128, 256x256, and 1024x1024 
grids and then propagated through turbulence. The width of the source was 
chosen so that the Fourier transform had approximately the same width on the 
final grid and the Gaussian sources on the different grids were normalized to 
the same energies. Since identical sources were applied in the spatial domain 
according to the analytical Gaussian formula, the spatial frequency 
representation of these sources improved as the grid size increased and the 
mesh became finer. The 1024x1024 grid served as an approximation to infinite 
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Figure 31 Estimates and least squares fits from spherical wave coherence 
length (triangles / asterisks) and HWHM (diamonds / pluses) of the maximum 
turbulence strength versus grid size for valid simulations. 



grid size. Due to its larger extent and finer mesh, the 1024x1024 spatial 
frequency representation of the Gaussians contained spectral energy in the 
region outside the spectral footprints of the smaller grids, as Fig. (32) illustrates. 
Since each grid started with the same total energy, the fraction of spectral 
energy in the 1024x1024 grid lying outside the spectral footprints of the smaller 
grids approximated the amount of energy aliased in the smaller grids. This 
process was then applied to propagations of the Gaussian beams through 
increasing strengths of turbulence, well into the significant aliasing regions for 



Figure 32 Illustration of the spectral representations of the Gaussian beam 
with different grid sizes. 
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the smaller grids. This method provided estimates of the fraction of energy 
aliased versus turbulence strength for the three smaller grids. 

To achieve the first goal of relating the amount of aliasing to a 
measurable simulation parameter, the average radial power spectral density of 
the propagated E-field was calculated for each turbulence strength with the 
smaller grids. Using these profiles, a spectral ratio was calculated, defined as 
the ratio of the power spectral density at the grid center to the power spectral 
density at the maximum grid radius. Figure (33) plots these results, and also 
includes spectral ratios for 512x512 and 1024x1024 grids. Using the 
information on fraction of energy aliased versus p 0 2 and spectral ratio versus 
p 0 2 , the turbulence strength p 0 2 served as the connection between fraction of 
energy aliased and the measurable simulation parameter of spectral ratio. 

Figure (34) shows plots of these spectral ratios versus corresponding fraction of 
energy aliased. The log-log plot behavior proved approximately linear over the 
range 0.001 to 0.1 , which appears very useful because 0.1% of energy aliased 
probably does not affect the simulation results significantly while more than 
10% of energy aliased probably does. 

Using the data points of Fig. (34) and performing a linear least squares 
fit to the logarithms 

log 10 (fl) = B log 10 (Ffl) ♦ log 10 4 (110) 
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or equivalently, 


( 111 ) 


= 10* {FR) B , 

where R refers to spectral ratio and FR refers to fraction of energy aliased, 
gives A= -1.1 ± 0.1, B= -2.1 ± 0.1 . Substituting the A and B values into Eq. 
(Ill) and rearranging, 


FR * 1 . = . 


( 112 ) 


This equation achieves the first goal of relating the fraction of energy aliased to 
a measurable quantity from the simulation. 

However, Eq, (112) should be used with caution because it appears to 
be specific to the narrow Gaussian source used to derive it. Figure (35) shows 
the same calculations carried out for a source that was the Fourier transform of 
an aperture of radius equal to 3/4 of the 64x64 grid radius (i.e. Airy-type 
source). The linear region did not appear, though the plot matched the 
Gaussian source plot when the fraction of energy aliased was greater than 
0.1 . The differences arise from the different spectral representations of the 
propagated E-fields. The abrupt, cylinder-shaped irradiance patterns start with 
more energy at higher spatial frequencies in the spectral domain than the 
Gaussian patterns and thus never get below 0.001 fraction of energy aliased 
even at very low turbulence strengths and high power spectral ratios. 






The second goal of determining the amount of aliasing that invalidates 
the simulation used the fraction of spectral energy aliased to generate another 
guideline for maximum strength of turbulence for a given grid size where the 
maximum fraction of energy allowed to be aliased in a simulation was 10%. 
Equation (112) predicts a spectral ratio of approximately 10 for 10% of energy 
aliased. Linearly interpolation between the data points of Fig. (33) provided 
estimates of the turbulence strengths C„ 2 (thus p 0 2 ) that gave a spectral ratio 
= 10, for 10% energy aliased, for the five grid sizes. Figure (36) plots these 
values and a least squares fit extended to larger grid sizes. The least squares 
fit as a function of grid size N was 

lOQioPo = 0.9 log 10 ( N) - 0.9, (113) 

or, 

= 0.1 A/ 09 . (114) 

The fraction of energy aliased also provided estimates of maximum p 0 2 
for the Airy-type source. Figure (37) plots the fraction of energy aliased versus 
turbulence strength for the 64x64, 128x128, and 256x256 grids. Again 
choosing 10% energy aliased and using linear interpolation yields data points 
for the solid line in Fig. (36) that shows the maximum turbulence strengths 
obtained from Fig. (37). For 10% energy aliased, the predictions for maximum 
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Figure 36 Least squares fit (pluses, dotted line) of maximum beta squared 
versus grid size for a Gaussian source that had 10% energy aliased and for 
Airy-type source (asterisks, solid line). 









P 0 2 for the three smaller grids using the Airy-type source agree within -20% with 
the predictions from the Gaussian source. Lower values of fractional energy 
aliased do not provide such agreement. 

Another measure of maximum turbulence strength involved the 
coherence length and grid size. Aliasing occurs because the E-field fluctuates 
significantly over scale sizes smaller than the grid element size. Intuitively, 
some connection should exist between coherence length of the E-field and the 
onset of significant aliasing. Fried's coherence length r 0 represents the distance 
over which the atmospheric MTF falls to exp(-3.44) = 0.032 (Fried, 1966, p. 
1380-1383). Equation (38), which relates p 0 2 to C„ 2 , and Eq. (66), which 
relates r 0 to C n 2 , and Eq. (92), which relates the grid element size Ax to N, 
when combined, yield the turbulence strength corresponding to a given grid size 
for a specific number (y) of Ax's per r 0 

„ , 0.677 y 5/3 A/ 5 ' 6 ,spherical wave 

PS = (US) 

1 0.629 y" 5 ' 3 A / 5 ' 6 ,plane wave 

Figure (38) plots these p 0 2 for y = 2, 2.5, and 3 for the spherical wave case. 

The 10% aliased energy line, indicated by plus symbols, lies nearest the 2.5 Ax 
per r 0 line. The E-field must be spatially sampled with 2.5 Ax per r 0 to limit the 
fraction of energy aliased to 10%, analogous to the Nyquist criterion of two 
samples per cycle. 
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Grid Size 


Figure 38 Theoretical spherical wave coherence length versus grid size 
requiring 2 (diamonds), 2.5 (triangles), or 3 (squares) grid elements per 
coherence length. 10% aliased energy line (pluses) also plotted for reference. 





Figure (39) provides a combined plot of the maximum turbulence 
strength p 0 2 that produced valid E-fields for a given grid size from all of the 
foregoing measures: intensity ratio, Strehl ratio, normalized irradiance variance, 
coherence length, HWHM, fraction of energy aliased, and coherence length 
(y = 2.5). Most estimates lie within a factor of 2 of each other, and the 10% 
energy aliased line (dotted/pluses) given by Eq. (114) represents an 
approximate lower bound to those estimates based upon onset of significant 
aliasing. 

Figure (39) provided three significant conclusions: (1) the computer 
simulations remained valid until approximately 10% of the energy became 
aliased (achieving the first goal), (2) using the 10% energy aliased line, the 
maximum turbulence strength p 0 2 for valid E-fields for a grid size N was 

Po^ = 0.1 A/ 09 , (116) 

and (3) maximum valid turbulence strength corresponded to approximately 
2.5 grid elements per r 0 (based on proximity of the y = 2.5 line to the 10% 
energy aliased line). Conclusion (2) implies that doubling the grid size N (with 
grid element size Ax given by Eq. (92)) slightly less than doubles the maximum 
P 0 2 that the grid can simulate, and that these investigations, which use a 
1024x1024 grid, should be valid up to P 0 2 msi> - 50. 
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ADDITIONAL SIMULATION PARAMETERS 


Additional aspects of a computer simulation were addressed to ensure 
validity, including the number of phase screens, the method of normalizing the 
irradiance variance, the number of realizations (propagations through 
turbulence) required for representative statistics, and the width of the final 
irradiance field on the grid. 

The number of phase screens must be large enough to represent the 
turbulence accurately along the path and produce proper irradiance and 
coherence statistics. Martin and Flatte (1988) determined the number of phase 
screens by requiring that the variance due to propagation over the distance Az 
between phase screens be less than 1/10 the variance from the total 
propagation over the distance L 

o?(A z) < 0.1 o*(L), (117) 

and additionally that the value of the variance from one step be less than 0.1 

of(A z) < 0.1 . (118) 

With these considerations, they generally used 20 phase screens for their 
simulations. 

These investigations examined the irradiance and coherence statistics 
directly and concluded that approximately 30 phase screens were required to 
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ensure simulation validity. Figure (40) shows a set of spherical wave 
normalized irradiance variance simulations for a 512x512 grid, p 0 2 = 1.5, and 
the number of phase screens varying between 2, 4, 8, 16, 32, and 64. 

Assuming the 64 phase screen case most closely approximated physical reality, 
as few as 8 phase screens gave a normalized irradiance variance within -10% of 
this value. By 32 phase screens, the normalized irradiance variance had 
stabilized to within 2% of the 64 phase screen value. For this level of 
turbulence that lies at the beginning of the saturation regime, the addition of 
phase screens beyond 32 did not affect the normalized irradiance variance 
significantly. 

Figure (41) shows the normalized irradiance variance from propagations 
with 2, 4, 8, 16, 32, and 64 phase screens at p 0 2 = 50 with a 1024x1024 grid 
and 30 realizations for each number of phase screens. This strength of 
turbulence represents the limit of simulation validity with the 1024x1024 grid 
and lies in the strong saturation regime beyond the normalized irradiance 
variance peak. In this case, 32 phase screens provided a normalized 
irradiance variance within 3% of the 64 phase screen case, and the 16 phase 
screen value was within approximately 10%. Because 32 phase screens 
appears to increase the accuracy of the simulation by -5% compared to 16 
phase screens, the larger number of phase screens was chosen for these 
investigations. 













Using the same 1024x1024 realizations, the coherence length r 0 and 
HWHM of the atmospheric MTF were calculated for the 2, 4, 8, 16, 32, and 64 
phase screen cases. Figure (42) plots the results. For coherence length, the 
32 phase screen value was within approximately 1% of the 64 phase screen 
value, and the 16 phase screen value was within approximately 3%, The 
HWHM plot indicates 2% and 3% agreements for 32 and 16 phase screens, 
respectively. The 32 phase screens used for these investigations thus proved 
sufficient for coherence length and HWHM calculations. 

The method of normalizing the irradiance variance and the number of 
realizations to use for statistical accuracy proved to be interrelated 
considerations. Turbulence diffracts and scatters energy outward from an 
initially well-defined beam. For the Airy-type source used in these 
investigations, the average irradiance over the central portion of the final 
propagated field was uniform for zero turbulence but became a function of 
radial distance from the propagation axis when turbulence was present. This 
was an artifact of the computer simulation that had to use a beam spatially 
confined to the grid rather than a true spherical (or plane) wave. Figure (43) 
shows irradiance versus radial distance r for the central 256x256 portion of 
1024x1024 grid simulations, averaged over 30 realizations, and using a 
turbulence strength p 0 2 = 50. The average radial irradiance varies by 
approximately 15% over the width of this calculation region. This radial 
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Figure 42 Coherence length (top) and HWHM (bottom) as a function of 
different numbers of phase screens at turbulence strength p 0 2 = 50 . 
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variation of the average irradiance was removed from the normalized irradiance 
variance calculations with a method similar to that of Flatte, Wang, and Martin 
(1993), The calculation region was restricted to half the radius of the zero 
turbulence irradiance pattern, which was previously chosen as one half the grid 
radius. With the 1024x1024 grid, this calculation region equaled the largest 
radius circle inside the central 256x256 portion of the grid. This disk was 
further divided into concentric rings and the irradiance in each ring averaged 
over all 30 realizations. These average ring irradiances were used to normalize 
the irradiance variance calculated from each field. Smaller ring size (thus more 
rings) reduced the number of points in the irradiance average for each ring and 
required more realizations to ensure a sufficient number of points to yield a 
stable average irradiance. 

The ring width had to be small enough to compensate for the radial 
variation in average irradiance and the number of realizations large enough to 
provide enough points to yield representative average values. To determine 
suitable ring width and number of realizations, 50 1024x1024 realizations with 
the Airy-type source at p 0 2 = 1.5, 3, 10, and 20 were run and the normalized 
irradiance variances calculated using ring widths of 128, 32, 8, 4, 2, and 1 grid 
elements, Ax, and with number of realizations in the average equal to 1, 5, 10, 
15, ... ,50. Figure (44) plots the results for the p 0 2 = 10 set. The pluses line 
corresponds to 1 ring of width 128 Ax (i.e. the whole calculation region) and a 
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stable normalized irradiance value was achieved using about 20 realizations. 
However, a single ring does not compensate for the radial variation in average 
irradiance and can thus give erroneous normalized irradiance values (for the 
same reason, it becomes sensitive to the peaking behavior from significant 
aliasing). The X line represents the division of the disk into 4 rings of width 32 
Ax and required ~25 realizations in the average to give a stable normalized 
irradiance variance. Smaller ring widths (lines with squares, diamonds, 
triangles, and asterisks lines) all showed similar behavior and required ~ 30 runs 
to reach a stable average normalized irradiance variance. The p 0 2 = 1.5, 3, 
and 20 runs all provided similar results. A similar p 0 2 = 10 series with a 
Gaussian source required a minimum of 8 rings (16 Ax each) and 30 
realizations to achieve stable average normalized irradiance variances. 
Consequently, these investigations used 30 realizations as a guideline for valid 
simulation, and chose a ring width of four grid elements to allow the greatest 
adjustment to real variations in average radial intensity and yet remain 
computationally efficient. As Fig. (44) shows, insufficient averaging over enough 
realizations to yield truly representative average values can reduce the 
normalized irradiance variance by 15 - 30%. 

Similarly, stable coherence length values (discussed in section, H. 
Coherence Length) required averaging over multiple realizations. To determine 
the number of realizations required in the average, 50 realizations using a 
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1024x1024 grid were generated for P 0 2 = 3 and 20 and the central 256x256 
portion of the fields again used to calculate the atmospheric MTFs. These 
individual MTFs were then averaged in groups of 1, 5, 10, 15, ... , 50 before 
calculating the coherence lengths. Figure (45) shows the resulting coherence 
lengths, r 0 . (These particular runs used L = 150 m and X = 420 nm, giving a 
Fresnel length = ja l / 2 it = 10 mm and resulting in millimeter-sized 
coherence lengths.) As few as 5 realizations in the average yielded coherence 
lengths within 4% of the 50-realization values. These investigations still used 
30 realizations in the average because these 30 fields were already available 
from the normalized irradiance variance calculations. 

The irradiance ratio provided an indicator of the appropriate final beam 
radius to use to avoid excessive scatter of energy off the grid. The maximum 
p 0 2 for a given grid varies with the final beam radius because larger radii scatter 
energy off the grid sooner and cause aliasing at lower p 0 2 . To characterize this 
behavior, 64x64, 128x128, and 256x256 grid simulations were run at turbulence 
strengths p 0 2 = [5x1c 4 , 5x10 2 ] for final beam radii of 4/8, 5/8, 6/8, and 7/8 grid 
radius, and the resulting irradiance ratios (average irradiance at center divided 
by average irradiance at grid radius) were plotted versus p 0 2 . Figure (46) 
shows the plot for the 256x256 grid. Again, the minima correspond to the onset 
of peaking and occur at lower p 0 2 for larger final beam radius. Linear least 
squares fit of the final beam radii to these 6 0 2 values for the minimum 
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irradiance ratios indicates that a final beam radius of approximately 0.7 grid 
radius corresponds with the 10% energy aliased cutoff p 0 2 for the grids. A final 
beam of 0.7 grid radius provided the largest illuminated central region while still 
meeting the 10% aliased energy criterion. To be somewhat conservative, these 
investigations used a final beam radius = 0.5 grid radius to reduce further the 
energy scattered off the grid. 

G. PHASE SCREENS 

1. Phase Screen Generation 

With the split-step method, the effects of turbulence along the 
optical path are introduced into the simulation by dividing the optical path into 
steps and applying a random phase to the complex E-field at each step. As 
long as the steps are small enough that geometrical optics approximately 
applies, the E-field only acquires a random phase change as it propagates 
across each step (Knepp, 1983). Diffraction as the field propagates across 
many steps then produces the amplitude variations. The random phases are 
assumed to be Gaussian distributed about a zero mean with variance 
proportional to the turbulence strength C n 2 and possessing spatial structure 
function consistent with the assumed Kolmogorov turbulence and appropriate 
inner scale. Knepp (1983) and Martin and Fiatte (1988) describe the process 
of generating the phase screen with these characteristics, as discussed below. 
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The phase screen generation begins in the spatial frequency 
domain by imposing the proper spatial structure function. An NxN grid of 
complex numbers ©o^.k,) is formed whose real and imaginary parts are each 
Gaussian distributed random numbers with zero mean and unity standard 
deviation. This represents the Fourier transform of a grid of 

uncorrelated Gaussian distributed random numbers 0 o (x p y) representing phases. 
The proper spatial structure function corresponding to turbulence statistics is 
imposed upon the random phases 0 o (x,y) by applying a filter A^.i^) to 

0 o (k*,s) 

eOc,,*,) = A(x x ,K y ) e 0 (K x , Ky ). (ii9) 

Taking the magnitude of both sides of Eq. (119), squaring, assuming that the 
filter function is real, and then taking expectation values gives 

( | 0( Kjr ,«C y ) | 2 > = A 2 (K x ,K y ) ( | 0 o (K x ,K y ) | 2 > = » 2 (K x ,K y ), (120) 

where use was made of the fact that the Gaussian random numbers 0 o (k x ,k v ) 
have a variance of 1. The two-dimensional power spectral density of the phase 
F s (K x ,Ky) is related to ©(k,,^) by (Goodman, 1985) 

F S ( K*.* y ) -< |e(K XI K y )| 2 )(Ax) 2 , (121) 

where Ak represents the grid element size in the spatial frequency domain. 

The Hankel transform of the power spectral density F s (k) gives the phase 
structure function D s (p) that characterizes the spatial distribution of the phase 
fluctuations of the E-field (Tatarski, 1961) 
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( 122 ) 


0 S (P) = /l 1 - ^.(KP) ] F s { X,0) K C/x, 
where local isotropy has been assumed. Tatarski derived the relation between 
the two-dimensional power spectral density of phase fluctuations F s (k,,k y ) and 
the three-dimensional spectrum of index of refraction fluctuations, <t>„(K), 

F s (x,0) = 2nk 2 L *„(x), (123) 

where k = y/x* + x y . This sequence of steps means that the phase screen 

can have the proper spatial statistics by starting with the proper spectrum of 
refractive index fluctuations (hence, the proper structure function). 

The spectrum of index of refraction fluctuations assuming 
Kolmogorov turbulence with inner scale is 

®„(k,z) = 0.033 C„ 2 (z) x’ 11 ' 3 F( kCq), ( 124 ) 

where F(k«„) gives the inner scale dependence for inner scale Cq (see Fig. 8). 
Substituting Eq. (124) into Eq. (123) gives the power spectral density of phase 
fluctuations 

F s (x,z) = Znk 2 L (0.033) C 2 (z) x 11/3 F (kCq), < 125 ) 

and using Eq. (121) specifies the proper form for ©(k^), the corresponding 
Fourier transform of the random phases 

< |e(K x ,K y ) I 2 ) = (Ax)- 2 2 nk z L (0.033) C 2 (z) x' 11 ' 3 F(xCq). ( 126 > 


110 



Equation (120) then gives the corresponding filter function to apply to 

the array of complex numbers 0 o (K„K y ) 

A(K x , Ky ) = (Ax)' 1 pnk 2 L (0.033) Cfo) k' 11/s F(kI 0 ). < 127 ) 

Since the Kolmogorov spectrum (oc k' 11/3 ) has a singularity at k = 0, the k= 0 
point in the filter function is set to zero, removing overall piston (i.e. common 
phase offset over whole screen) from the phase (Cochrane, 1985) and keeping 
the spectral energy finite. Conversion to a discrete grid representation occurs 
by substituting = n x Ak, = n y Ak, k 2 = (Ak ) 2 (n x 2 + n y 2 ), and 
Ak = 27 t/(N Ax), where Ax is the grid element size in the x-y domain and (n„,n y ) 
are grid coordinates. The Fourier transform (FT) of the filtered array of random 
variables gives the phase screen in the spatial domain 

0 (x,y) = 0.0984 k sJC*{z) L (N Ax) 5 ' 8 FT[ 9 0 (n x ,n y ) ]. < 128 ) 

Since this phase screen is actually complex-valued, both the real and imaginary 
parts represent valid random phase screens that were separately applied to the 
E-field (Cochrane, 1985). 

2. Low Spatial Frequency Correction 

Due to the finite size of the grid, the above phase screen will not 
have the proper structure function for separations of the order of the grid width. 
Low spatial frequency components, especially tilt-type terms, are under- 
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represented (Cochrane, 1985). These computer simulations incorporate an 
algorithm formulated by Cochrane that employs an expansion of the phase 
screen in a Karhunen-Loeve basis set (whose components are mutually 
orthogonal) to correct some of the low spatial frequency terms. 

The general idea of the low frequency correction builds upon the 
above procedure to generate a phase screen. The low spatial frequency 
contribution to the phase screen is a superposition of orthogonal low spatial 
frequency terms, just like the superposition of complex exponential terms by 
Fourier transform in the phase screen generation process above. The strength 
of each low spatial frequency term is a Gaussian random variable with an 
appropriate variance, just as the Gaussian random numbers above were filtered 
to give the proper variance and thus determine the strength of the 
corresponding complex exponential in the spatial domain. Thus, the two 
objectives involve finding an appropriate set of orthogonal low spatial frequency 
functions, and determining the corresponding variances applicable to 
atmospheric turbulence. 

An arbitrary function can be expanded in terms of Karhunen- 
Loeve functions that are orthogonal by definition. To determine a set of 
Karhunen-Loeve functions appropriate for Kolmogorov turbulence, Cochrane 
builds on the work of Noll and (1975) considers expansion of an arbitrary 
function <t>(r,0) over a circular aperture of radius R in terms of Zernike 
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polynomials (Born and Wolfe, 1970) 


<!>(/■, 0) = 'Za ] Z J { p,e), (129) 

/= 1 

where p represents the normalized distance r/R, a j are the expansion 
coefficients, and Z t are the Zernike polynomials. With W(r/R) representing the 
aperture function, the coefficients are 

aj = (1 IR*) jd 2 r W(r/R) 4>(r,0) Zj{rlR,Q). (130) 

Noll assumed that these coefficients were Gaussian random variables with zero 
mean and with a covariance 

( a;a y ) = [dp [dp ' W(p) W) ZJL P.0) < <Kflp) Wp') > Z/p'.e'j.OSl) 

Fourier transforming to the spatial frequency domain 

(a/ay) = ff di. d% Qj(%) ® (k/R, k'/R) Qy{*') , (132) 

where Qj(k) represents the Fourier transform of the jth Zernike polynomial, and 
0 (k7R,k 7R) represents the Kolmogorov spectrum of phase fluctuations. Noll 
analytically performed the integrals to give a covariance matrix in which the 
terms represent the expected covariances due to Kolmogorov turbulence. 

Cochrane (1985) notes that the Zernike polynomials cannot be 
used to form an orthogonal expansion of the turbulence-distorted phase 
because the expansion coefficients are correlated, indicated by nonzero off- 


113 



diagonal elements in Noll's covariance matrix. However, the eigenvectors of 
the Zernike covariance matrix serve as a Karhunen-Loeve basis set. These 
eigenvectors K p can represent turbulence because they are not correlated, i.e. 
each eigenvector K p is formed by superposition of Zernike polynomials in such 
a way that the K p are orthogonal, satisfying the first objective. Additionally, the 
corresponding eigenvalue X p multiplied by (D/r 0 ) 5/3 (where D is the aperture 
diameter and r D is Fried's coherence length (Fried, 1965)) gives the appropriate 
variance for that K p spatial component corresponding with Kolmogorov 
turbulence, satisfying the second objective. 

Specifically, the low spatial frequency contribution to the phase 
screen can be expanded in terms of these Karhunen-Loeve components K p 
(Cochrane, 1985) 

Pmm 

4>(r)= E Y,K,(-^), (133) 

p =1 

where p ma)t represents the number of low spatial frequency Karhunen-Loeve 
terms included and the coefficients y p are Gaussian random numbers with 
variance X p and scaled by (D/r 0 ) 5 ' 3 to the specific strength of turbulence used. 
The simulations use the first five terms ( p maK = 5 ) as a compromise between 
completeness of low spatial frequency correction and computational efficiency. 
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Figure (47) shows a surface plot of one realization of the first two terms of the 
correction, which are very close to x- and y-tilt (i.e. phase terms linear in x and 
y, respectively). Figure (48) shows a realization of the 3rd, 4th, and 5th 
correction terms, which resemble the wavefront aberrations associated with 
defocussing and astigmatism in a conventional imaging system. To implement 
the low spatial frequency phase correction, (1) the scalar product was formed 
between the initial Fourier transform phase screen and each Karhunen-Loeve 
function K pt giving the relative strength of that K p in the initial phase screen; (2) 
this amount of each spatial component Kp was then subtracted from the phase 
screen; and (3) the K„ component was then added back to the phase screen in 
the proper amount given by the product of a Gaussian random number y p with 
variance k p and the factor (DlrJ 513 to scale to the particular strength of 
turbulence used. 

Cochrane's computer routines only calculate the Karhunen-Loeve 
correction terms over the largest circle that fits inside the calculation grid, as 
shown in Figs. (47) and (48). Correction of the E-field over the entire computer 
simulation grid requires the Karhunen-Loeve terms be calculated over an area 
that is at least JZ larger on each side than the field grid (Cochrane, 1985). If 
the E-field grid size N is chosen as a power of 2, then the Karhunen-Loeve grid 
must be 2Nx2N. These simulations however used only an NxN grid 
(1024x1024) because the Sun SparcStations could not handle the memory 
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Figure 47 Terms 1 and 2 of Karhunen-Loeve low spatial frequency correction 
to phase screen, represented in optical path length for wavelength = 500 nm. 
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Optical path length (m) 



Figure 48 Terms 3, 4, 5 of Karhunen-Loeve low spatial frequency correction to 
phase screen, represented in optical path length for wavelength = 500 nm. 
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requirements of the doubled grid (2048x2048) without a major revision of the 
propagation code. Since the source function was chosen to minimize energy 
scattered off the grid, very little energy fell in the uncorrected corners of the grid 
and the simulation remained valid. 

Noll's covariance matrix assumes Kolmogorov turbulence 
spectrum with zero inner scale making the Karhunen-Loeve correction terms 
strictly apply only to this spectrum. These simulations however incorporate 
nonzero inner scales into the spectrum of refractive index fluctuations. 

Because the first five Karhunen-Loeve terms cover scale sizes on the order of 
the grid width, which is much larger than the inner scale size the Karhunen- 
Loeve correction terms derived for zero inner scale remain valid for correcting 
nonzero inner scale phase screens. 

Cochrane (1985) showed that such low spatial frequency 
correction greatly improved the phase structure function. Two terms corrected 
the structure function to within 10% of the theoretical Kolmogorov behavior, and 
five terms corrected to within 5%, compared with 30 - 1000% discrepancies 
without any correction at low spatial frequencies. Figure (49) plots the 
normalized irradiance variances from computer simulations as a function of p 0 2 
for zero Karhunen-Loeve low frequency correction, two-term correction, and 
five-term correction. Zero correction underestimated the irradiance variance in 
the Rytov regime by approximately 5%, and the two-term tilt-type correction also 
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Figure 49 Normalized irradiance variance, normalized by p 0 2 , versus 
turbulence strength for zero (dotted), two-term (dashed), and five-term (solid) 
Karhunen-Loeve low spatial frequency corrections to phase screens. 
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underestimated by about 5%. The five-term correction with its focussing type 
terms raised the variance to within about 2% of the theoretical variance. In the 
saturation regime, the two- and five-term corrections fit better. Figure (50) plots 
the coherence length r 0 , normalized by the theoretical coherence length, as a 
function of turbulence strength for zero, two-, and five-term Karhunen-Loeve low 
spatial frequency corrections. In the Rytov regime where simulation should 
closely approximate theory, the zero correction overestimated the coherence 
length by approximately 35%. The tilt-type correction (two terms) estimated the 
coherence length within about 5%, and the five-term correction achieved 
agreement within about 2%. These behaviors formed the guideline that some 
type of low spatial frequency correction was required to achieve valid 
coherence lengths from computer simulation. 

This low spatial frequency correction method remained 
computationally feasible because the Noll covariance matrix only needed to be 
calculated once and because only a few terms of the Karhunen-Loeve 
expansion were used. However, the code can become memory intensive 
because it saves a full NxN grid for each Karhunen-Loeve function in addition 
to the NxN phase screen itself. Implementation with 2048x2048 or larger grids 
becomes problematical except on very large computers with about one gigabyte 
of RAM. Fried(1993) has proposed a simpler x- and y-tilt correction algorithm 
and indicates that this performs almost as well as including the higher order 
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Figure SO Spherical wave coherence length, normalized by theoretical 
coherence length, versus turbulence strength for zero (dotted), two-term 
(dashed), and five-term (solid) Karhunen-Loeve low spatial frequency 
corrections to phase screens. 
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Karhunen-Loeve terms, while being less memory intensive and faster than the 
Karhunen-Loeve correction method. Figs. (49) and (50) indicate that the tilt- 
only correction may systematically underestimate the normalized irradiance 
variance by about 5% and overestimate the coherence length of a spherically 
diverging wave by about 5% in the Rytov regime. Depending on the 
application, the Karhunen-Loeve correction to higher order terms may prove a 
useful refinement. 

H. COHERENCE LENGTH 

The coherence length of the E-field was calculated from the atmospheric 
MTF with the theory outlined in Chapter II, but the actual implementation of the 
MTF calculation and parameterization of the coherence length r 0 required 
careful consideration. The calculation methods that worked best for these 
investigations are discussed in this section. 

Equation (60) of Chapter II relates the atmospheric MTF and the 
spherical wave structure function to the coherence properties of the E-field 

MTF ■ = e' 1 * m) = W-v-p-v = (< U(n UV't) > > (i 34 ) 

( W(F) W(F+f) ) ' 

Again, U(?) represents the E-field and W(/) represents the aperture function. 
The autocorrelations of U and W indicated in Eq. (134) could be done by 
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summing the complex products of the E-fields over multiple pairs of points or by 
implementing the autocorrelation via FFT techniques. Both versions were tried 
for these investigations and produced similar results, but the FFT version 
provided a much more thorough autocorrelation with greater computational 
efficiency. 

The FFT autocorrelation technique used the MCF introduced earlier. 
Equation (57) defined the mutual coherence function as 

MCF rn < U*(t„t,) U(t 2 ,Q ), (135) 

which, for a single time t, = t 2 = t and assuming homogeneity, may be written 
as the spatial autocorrelation of the E-field 

MCF{t") = f dt U’{!) U(t + t"). (136) 

Substituting Eq. (136) into the Fourier transform identity, 

MCF{t>) = / df[ [ dt" MCF(t") er' 2 *''"] (137) 

gives 

MCF(F) = j dt \ f dt"[ f Cf*(/) U(M") ] e~ ,2 ' t - r " ] e*' 2 *'''. (138) 

Rearranging the integrations, 

MCF(f) = f df [ / dt U'(f) [ j dt" U(t+t") (139) 

Changing variables to s = t + t", t" = $ - ?, dt" = dS, 
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MCF(F) = ( df [ / df U'(f) [ [ dS U($) e-' z ' ,a ] e * l2 ' tf ] e^ hf '. ( 14 °) 

The inner integral is the Fourier transform of U(l) and is denoted by FT[ U ]. 

The next innermost integral is the inverse Fourier transform of U*(?) and is 
denoted by IFT[ U* ]. Then, 

MCF(f) = / df FT[ U] IFT[ L/* ] e^ !t ', (W) 

which is yet another inverse Fourier transform, symbolically written 

MCF(f) = IFT[ FT[U] IFT[U'] ). ( 142 ) 

Equation (142) expresses the autocorrelation (MCF) of the E-field in terms of 
Fourier transform techniques easily implemented with discrete Fourier 
transforms. 

This Fourier transform method of calculating the autocorrelation of a 
function is faster and more complete than the more laborious technique of 
averaging the products of the E-field at point pairs. The product E*( )E( t 2 ) is 

complex, but because of the averaging that occurs in the autocorrelation, the 
real part attains a stable value while the imaginary part averages toward zero. 

For the FFT implementation on a 1024x1024 grid, the imaginary part is typically 
~10' 9 while the real part ranges between 0 and numbers on the order of unity. 

The point pairs technique produces an equivalent real part but reduces the 
imaginary part down to only -10' 3 . These investigations used the FFT version for 
all autocorrelations. Additionally, the autocorrelations were normalized by their 
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zero lag values to ensure that the magnitude of the MTF atmos calculated by 
Eq. (134) lies between 0 and 1. 

Because the spherical wave structure function refers to points across a 
spherical wavefront, U and W must similarly represent the E-field across a 
spherical surface. However, the computer simulation implements the E-field on 
a plane perpendicular to the axis of propagation and incorporates the spherical 
wave nature of the E-field by applying a quadratic phase curvature across the 
plane. To convert from this plane representation of the E-field in the simulation 
to a spherical surface representation appropriate for the autocorrelation 
calculations of Eq. (134), the quadratic phase factor across the plane must be 
removed from the E-field. This effectively assumes that the amplitude and the 
phase fluctuations of the E-field across the plane of the computational grid 
closely approximate the E-field across the spherical wave. This assumption is 
justified because the maximum physical separation of the true spherical 
reference surface from the plane surface of the grid is at most 1/50 the 
coherence length of the E-field. However, the removal of the quadratic phase 
curvature from the E-field of the grid proves crucial to the autocorrelation of U 
because the quadratic phase factor undergoes ~130 multiples of 2n phase 
change between the center and the outer edge of the grid for these simulations 
and would otherwise completely obscure the actual E-field fluctuations. 
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Propagation with a plane or beam wave requires a different approach to 
remove the proper amount of phase curvature. For a plane wave, the 
wavefront coincides with the plane of the grid so that no phase curvature needs 
to be removed. A pure beam wave (Gaussian profile) exhibits a spherical 
phase with a radius of curvature larger than the propagation distance L, and 
this phase could be calculated analytically and removed. However, the use of a 
finite E-field confined to the propagation grid introduces phase effects other 
than simple quadratic curvature. Fortunately, the E-field propagated through 
zero turbulence contains this phase curvature information (Walters, 1994). For 
spherical wave propagations, the removal of the phase curvature by the 
analytical calculation and by using the zero turbulence propagated field 
provided identical coherence lengths. These investigations used the latter 
method for the coherence length calculations for all beam-like and plane-wave- 
like propagations. 

Equation (134) requires that the autocorrelation of U must be averaged 
over multiple realizations to achieve the long exposure MTF. To implement this 
requirement, the autocorrelations from several realizations were calculated and 
averaged together, and a coherence length then determined via Eq. (134). This 
method produced coherence lengths that agreed with theory within ~5% in the 
Rytov regime. Recalling Fig. (45) for the coherence length versus number of 
realizations used in the average, the number of fields included in the MTF 
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average may be as few as 5, though 20 provided a more statistically 
reproducible coherence length. Again, these simulations used 30 realizations 
because these fields had already been generated for the normalized irradiance 
variance calculations. 

Figure (51) plots the coherence lengths calculated from the average MTF 
and also plots the average of coherence lengths calculated with single 
realization MTF's. The average of individual realization coherence lengths 
exceeded the averaged MTF coherence length by -20% at low turbulence 
strengths but eventually agreed within 5% near the saturation regime. The 
single realization MTF's showed a larger coherence length at low turbulence 
strengths because the contributions from low spatial frequency components had 
not been reduced through averaging. At these low turbulence strengths, 
multiple realizations were required to average these low spatial frequency 
contributions and to achieve the appropriate long exposure MTF. At higher 
turbulence strengths near saturation, the E-field had more energy at high spatial 
frequencies that dominated the MTF. The autocorrelation for a single 
realization now averaged over many coherence lengths and yielded an MTF 
close to the average MTF for multiple realizations. 

Using a finite beam to approximate a spherically diverging wave 
produced a radial dependence of the average irradiance that affected the 
coherence length calculation, but only to a small degree. To compensate for 
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Figure 51 Coherence lengths r 0 obtained by averaging r 0 's from individual 
realizations (diamonds) and from averaging the MCF's of individual realizations 
prior to calculating r 0 (triangles). 
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this radial dependence, each E-field was divided by the radial average E-field 
magnitude from the 30 realizations. Similar to the irradiance variance 
calculation, this average E-field was calculated by dividing the grid into rings 
one grid element wide, taking the square root of the average intensity for each 
ring over the 30 realizations, and then performing an area weighted, running 
mean across the rings to smooth the variations. While this radial compensation 
proved essential for the normalized irradiance variance calculation, it only 
changed the coherence length by -1%, which was significantly less than the ~5% 
discrepancy from the low spatial frequency correction. 

The above considerations allowed calculation of the right-hand side of 
Eq. (134) for the atmospheric MTF. Because of the statistical nature of the 
propagation through turbulence, no set of realizations yielded an atmospheric 
MTF that exactly followed the exponential rolloff with distance predicted by the 
left-hand side of Eq. (134). Methods had to be developed to parameterize 
these atmospheric MTF's from the simulations and extract an appropriate 
coherence length corresponding to the structure function D(r). 

For the case of Kolmogorov turbulence, Fried (1966, p. 1380-1383) 
derived the wave structure function and expressed it in terms of the single 
coherence length parameter r 0 



where n = 5/3. The atmospheric MTF then becomes 



To extract the coherence length r 0 and the exponent n (allowing the possibility 
that n may vary), take the natural logarithm of both sides and rearrange 
(Walters, 1993) 

■ds ,n( " rF -“ M) '(•£)"• (145) 

Taking the natural logarithm again, 

ln( ‘ 3^4 " '"M " n ln ( / o)- < 146 > 

This has the linear form y = ax + b, where y represents the left hand side, 
a = n, x = ln(r), and b = -n ln(r 0 ). The atmospheric MTF can now be 
characterized with two parameters, r 0 and n, or just the single parameter r 0 
assuming n = 5/3. 

To implement these parameterizations, the atmospheric MTF was first 
calculated from the E-field using the autocorrelation methods above and then 
radially averaged to yield MTF almoa (r) for use in Eq. (146). A linear least 
squares fit calculation provided the slope a = n and the intercept n ln(r 0 ), giving 
r 0 . To obtain the single parameter characterization, n was set to 5/3 in 
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Eq. (146) and least squares techniques applied to obtain the intercept and thus 
the single parameter r 0 . 

Figure (52) shows the coherence lengths calculated with these least 
squares methods, divided by the Rytov theory coherence length for Kolmogorov 
turbulence. The two parameter characterization with both exponent n and 
coherence length r 0 never provided a consistent, smoothly varying coherence 
length. At low turbulence strengths in the Rytov regime where the E-field 
coherence length is larger than the calculation aperture, the two parameter fit 
predicted coherence lengths up to 50% higher than the theoretical value and 
thus appears unreliable. Additionally, it showed an anomalous bump around 
(3 0 2 ~ 0.5 . The single parameter least squares technique with n = 5/3 
accurately characterized the coherence length within 5% at low turbulence 
strength, but still showed the bump at p 0 2 - 0.5 . 

An alternate technique to obtain a single parameter characterization 
assumed n = 5/3 in Eq. (144) and used a binary-type search, or iterative fit, to 
find the coherence length r 0 that minimized the variance between the average 
MTF atm05 (r) and the right-hand side of Eq. (144). Though not analytical, the 
resulting r 0 gave an MTF that often fit the actual MTF atmoa (r) more closely by eye 
than the least squares methods, especially for low turbulence where coherence 
lengths were larger than the grid size. To implement the technique, an initial 
large range of r 0 (for example, 0 to 100 m) was divided in half, the midpoint r 0 's 
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Figure 52 Coherence lengths calculated by one-parameter (dashed line) and 
two-parameter (dotted line) least squares techniques, the iterative fit technique 
(solid line), and the HWHM (dash-dot line). 




of each half were substituted into Eq. (144), and the variances were calculated. 
Whichever r 0 provided the least variance became the middle of the next range, 
and the other r 0 became the new high or low boundary. The process was then 
repeated, so that each iteration reduced the range of possible r 0 by 1/3. This 
procedure was iterated 100 times or until the variance was less than 1x1 O' 7 . 

Figure (52) shows the coherence lengths predicted by this iterative fit 
method along with the coherence lengths from the one- and two-parameter 
least squares fits, all divided by the Rytov theory coherence length for 
Kolmogorov turbulence. The iterative fit coherence length was -5% too large at 
low turbulence but did not show the anomalous bump around p 0 2 - 0,5 . 

In the saturation regime where turbulence is high and/or path lengths are 
long, multiple scattering becomes significant and saturates the irradiance 
variance (Martin and Flatte, 1988). Correspondingly, this physical phenomenon 
may also affect the coherence length and modify the value of n or the form of 
Eq. (144) for the saturation regime. The half width at half maximum (HWHM) 
provided a coarser, one parameter method of characterizing the atmospheric 
MTF that did not depend on any assumptions about the form of the structure 
function and could be calculated directly from the atmospheric MTF by linearly 
interpolating between the pair of points bounding MTF atmo ,(r) = 1/2. Figure (52) 
plots the HWHM, divided by the HWHM predicted by Rytov theory assuming 
Kolmogorov turbulence. The HWHM plot starts at p 0 2 = 0.05 because lower 
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turbulence strengths gave a coherence large enough that the MTF atmos (r) had 
not reached half its maximum value within the calculation region. The HWHM 
lengths followed the theoretical values within 5% in the Rytov regime and did 
not show the anomalous bump around p 0 2 - 0.5 . Note that the iterative fit 
coherence lengths agreed with the HWHM plot better than the least squares 
techniques. The HWHM and iterative fit r 0 proved to be the most stable 
parameterizations of the E-field coherence length, and were used in all 
subsequent coherence length comparisons and plots. 

Other factors were also considered in the coherence length calculation. 

As stated earlier, some choices for point source, such as the Airy-type source, 
required more energy at high spatial frequencies than others, such as a narrow 
Gaussian. While this difference produced < 2% effect in the normalized 
irradiance variance calculations, it appeared to have more effect on coherence 
length calculations. Figure (53) shows that an Airy-type source produced 
coherence lengths up to 10% higher at low turbulence strengths than those 
from the Gaussian-type source of Eq. (109). Therefore, only the latter type 
source was used for further investigations. 

For an arbitrary spectrum of refractive index fluctuations <t> n (ic,z), the 
integral formulation Eq. (61) provided the wave structure function for the 
atmospheric MTF in Eq. (134). Specifically, numerical integration of Eq. (61) for 
the spectra with grid cutoff, Gaussian, and Hill/Frehlich viscous convective 
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Figure 53 Iterative fit coherence lengths versus turbulence strength for Airy- 
type source (dotted), Gaussian-type source (solid), and Gaussian-type source 
with zero Karhunen-Loeve low spatial frequency corrections (dashed). 
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enhancement inner scales allowed comparison of theory (Eq. (134)) and the 
simulation MTF almos (r) (described in Chapter IV). However, in the structure 
function numerical integration, care had to be exercised in properly treating the 
low spatial frequency portion of the integrand since the spectra contained a 
k' 11 ' 3 singularity as k approached zero. This low frequency portion was critical 
to obtain coherence lengths that approached the Kolmogorov theoretical values 
in the Rytov regime. The integral was successfully evaluated by integrating 
analytically for 0 < k < K min , ( where x mjn ~ IxIO^m' 1 ) (Walters, 1994) since the 
inner scale function F(k{ 0 ) = 1 here and this portion of the spectrum remains 
Kolmogorov. The remaining portion of the integral that contained the inner 
scale contribution was carried out numerically. More realistic spectra could also 
have included an outer scale, but again no universal form of outer scale exists 
due to anisotropy of the atmosphere at large scale sizes. When included in the 
numerical integrations for test purposes, an outer scale raised the coherence 
length compared to the Kolmogorov case. However, these investigations did 
not use an explicit outer scale in the simulations. 
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IV. RESULTS 


A. NORMALIZED IRRADIANCE VARIANCE 

The computer simulation guidelines and considerations discussed in 
Chapter III were implemented to investigate the behavior of the normalized 
irradiance variance and E-field coherence length in the Rytov and saturation 
regimes for the grid cutoff, Gaussian, and Hill/Frehlich viscous-convective 
enhancement inner scales. Specifically, the simulations apply to stratospheric 
propagation with propagation distance L = 200 km, wavelength X = 500 nm, 
strengths of turbulence p 0 2 = [5x10 J \ 50], and inner scale sizes [0, 15] cm. 

The simulations used a 1024x1024 grid with grid element size given by Eq. 

(92), an Airy-type source modified to produce a final zero turbulence irradiance 
pattern with edges apodized by a Gaussian (Eq. (109)) and width corresponding 
to half the grid width, 32 phase screens utilizing a five-term Karhunen-Loeve 
low spatial frequency correction, and 30 realizations in each set of runs. The 
central 256x256 portion of each propagated E-field was used for the normalized 
irradiance variance and coherence length calculations. 

For the Gaussian inner scale values of 0 (grid cutoff), 5, 10, and 15 cm, 
Fig. (54) plots the normalized irradiance variance versus turbulence strength p 0 2 
in the Rytov regime from both numerical integration of the equation from Rytov- 
Tatarski theory, Eq. (34) (dotted lines), and from computer simulations 
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Figure 54 Normalized irradiance variance (divided by (3 0 2 ) from simulation 
(solid) and numerical integration (dotted) for Gaussian inner scales of 0, 5, 10, 






(solid lines). All values are normalized by p/ . Numerical integration values 
for 0 cm inner scale agreed within 1% of the theoretical zero inner scale values. 
The difference resulted because the numerical integration was limited to spatial 
frequencies below the grid cutoff K ma> . Larger grid cutoff values gave closer 
agreement. The simulation normalized irradiance variances agreed within 2% 
of the numerical integration values for all four inner scale values examined. 
Nonzero Gaussian inner scales reduced the normalized irradiance variance 
below the zero inner scale value (by 10%, 25%, and 40% for the 5, 10, 15 cm 
cases, respectively). Intuitively, the finite inner scale suppressed the higher 
spatial frequency index of refraction fluctuations and thus reduced the variance. 

Figure (55) plots the normalized irradiance variance (divided by p 0 2 ) for 
the single turbulence strength p 0 2 = 5x10" 4 and Gaussian inner scale sizes of 0 
(grid cutoff), 5, 10, and 15 cm. The numerical integration values (dotted line) 
and computer simulation values (solid line) showed an almost linear decrease 
of the normalized irradiance variance with increasing inner scale size in the 
Rytov regime. As Flatte, Wang, and Martin (1993) point out, the Gaussian 
inner scale does not accurately describe the inner scale observed in the 
atmosphere but retains usefulness because it facilitates some theoretical 
calculations. 

For the Hill viscous-convective enhancement inner scale sizes of 0 (grid 
cutoff), 5, 10, and 15 cm, Fig. (56) plots the normalized irradiance variance in 
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Figure 55 Normalized irradiance variance from simulation (solid) and numerical 
integration (dotted) for Gaussian inner scales of 0,5,10,15 cm at p 0 2 = 5x1 O' 4 . 
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Figure 56 Normalized irradiance variance (divided by p 0 2 ) from simulation 
(solid) and numerical integration (dotted) for Hill viscous-convective 
enhancement inner scales of 0,5,10,15 cm. 




the Rytov regime versus turbulence strength p 0 2 . Numerical integration of the 
Rytov-Tatarski theory, Eq. (34) (dotted lines), and computer simulation (solid 
lines) agreed within 2%. For the smaller inner scales, the normalized irradiance 
variance exceeded the zero inner scale values (by 30% for the 5 cm case). 

The viscous-convective enhancement of the strength of higher spatial 
wavenumber fluctuations near the inner scale wavenumber increased the 
variance. Yet, for large enough inner scale, the rolloff beyond the enhancement 
suppressed the higher spatial frequency fluctuations enough to eventually 
reduce the variance below the zero inner scale variance (by 30% for the 15 cm 
case). The 10 cm values happened to lie within 1% of the 0 cm values. 

Figure (57) plots the normalized irradiance variance (divided by P 0 2 ) for 
the single turbulence strength p 0 2 = 5x10"* and the Hill viscous-convective 
enhancement inner scale sizes of 0 (grid cutoff), 2, 3, 4, 5, 6, 7, 10, and 15 cm. 
Numerical integration of the Rytov-Tatarski results (dotted line) and computer 
simulation (solid line) clearly illustrated the rising and then falling behavior of 
the normalized irradiance variance with increasing inner scale size in the Rytov 
regime. The normalized irradiance variance achieved a maximum for ^ ~ 4 cm, 
which was about 30% of the Fresnel length R, = -jX L / 2 n = 12.6 cm. 

The dashed line in Fig. (57) shows simulation values using the Frehlich 
parameterization of the viscous-convective enhancement inner scale. The 
Frehlich inner scale shifted the plot slightly to smaller inner scale sizes, 
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maximized at 2% less than the Hill maximum, and matched the Hill values 
within 3% over the range of inner scale plotted. Additionally, simulation runs 
using 4 cm Hill and Frehlich inner scales agreed within 3% over the range of 
turbulence strengths 5x1 O' 4 < p 0 2 < 50. Thus, the Hill and Frehlich versions 
of the viscous-convective enhancement inner scale perform almost identically. 

Previous investigations have illustrated the dramatic monotonic rise in 
normalized irradiance variance in the saturation regime as the inner scale size 
increases (Martin and Flatte, 1988). Figure (58) shows normalized irradiance 
variance from computer simulations for 0 (grid cutoff), 5, 10, and 15 cm 
Gaussian inner scales and a propagation path of 200 km. Figure (59) shows a 
corresponding plot for 0 (grid cutoff), 5, 10, and 15 cm Hill inner scales. The 
turbulence values range from the Rytov regime (low turbulence with p 0 J < 1) to 
the saturation regime (high turbulence, 3 0 J £ 1 ). The Rytov regime showed 
again the behaviors illustrated with Figs. (54) - (57). Increasing Gaussian inner 
scale size produced monotonically decreasing normalized irradiance variance. 
Martin and Flatte (1988, 1990) provided similar plots of normalized irradiance 
variance with a Gaussian inner scale. The Hill viscous-convective 
enhancement caused the normalized irradiance variance to rise and then fall as 
the inner scale size increased. However, in the saturation regime, the 
normalized irradiance variance increased monotonically with increasing inner 
scale size for both the Gaussian and Hill inner scales. The transition in 
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Figure 58 Normalized irradiance variance over Rytov and saturation regimes 
for Gaussian inner scales of 0 (solid), 5 (dashed), 10 (dash-dot), and 15 cm 








behavior occurred with the onset of saturation around p 0 2 ~ 1 . This crossing 
behavior has been plotted for the log intensity variance with the viscous- 
convective inner scale by Hill and Clifford (1978). 

Figures (54) - (57) illustrate the close agreement between numerical 
integration and computer simulation values for the normalized irradiance 
variance at low strengths of turbulence. This agreement provided a validity 
check on these computer simulations that incorporated an inner scale. 

B. COHERENCE LENGTH 

Coherence lengths of the E-field were calculated from the same 
1024x1024 simulation runs used for the normalized irradiance variance 
calculation. The average atmospheric MTF was formed from 30 realizations 
using the FFT autocorrelation method, and then the corresponding coherence 
length r 0 and HWHM were calculated. The iterative fit r 0 's were used for 
comparisons because they most closely followed the HWHM behavior. The 
HWHM may provide the coarsest measure of coherence length but, since it 
requires no assumptions about the form of the MTF, it may also be the most 
reliable. 

Figure (60) shows the coherence length r 0 from numerical integration of 
Eq. (61) for an approximately zero inner scale and normalized by the theoretical 
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Figure 60 Coherence length from numerical integration for zero inner scale 
versus turbulence strength. Values normalized by Kolmogorov turbulence ze 
inner scale coherence length. 



Kolmogorov turbulence zero inner scale coherence length r 0 of Eq. (66). This 
plot indicates that the numerical integration coherence lengths calculated in 
these investigations were accurate within 2%. 

Figure (61) shows the coherence length r 0 from numerical integration of 
Eq. (61) with the grid cutoff K ma)( = 318 rad/m for the 1024x1024 grid. The grid 
cutoff did not affect the coherence length until p 0 2 « 1.5. Above that level of 
turbulence, the absence of the very high spatial frequency contribution to the 
integral caused the coherence length to become larger than the theoretical zero 
inner scale value. The grid cutoff inner scale implicit in the computer 
simulations with the finite grid should have a similar effect at these high 
turbulence values. 

Figure (62) shows the coherence length r 0 , normalized by the theoretical 
coherence length, from computer simulation of a spherically diverging E-field. 

In the Rytov regime, the simulation agreed with theory to within the ~5% 
overestimation due to using only five terms in the Karhunen-Loeve low spatial 
frequency correction to the phase screens (see Fig. (50)). In the saturation 
regime, the simulation coherence length (solid line) dropped below the 
theoretical prediction (by -25% at p 0 2 = 50) and the HWHM (dotted line) mirrored 
this decrease. For strong turbulence, the first order perturbation theory basis 
for coherence length appears to lose validity, as it did for the normalized 
irradiance variance in the saturation regime. 
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Figure 61 Coherence length from numerical integration for grid cutoff inner 
scale versus turbulence strength. Values normalized by Kolmogorov turbulence 
zero inner scale coherence length. 
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Figure 62 Coherence length r 0 (solid line) and HWHM (dotted line) for 
spherical wave propagation through turbulence (Values normalized by 
theoretical coherence lengths). 





The drop in coherence length presumably came from strong scattering in 
high turbulence, but a possible origin in the computer simulation was also 
considered. The simulations started with an approximate point source to 
emulate spherical divergence of the E-field and required the E-field to remain 
basically confined within the simulation grid over the propagation. The resulting 
E-field properties could have differed from those of a true spherical wave. To 
investigate this possibility, the width of the final irradiance pattern was varied 
such that wider final irradiance fields (hence narrower sources) more closely 
approximated a true point source. Figure (63) plots the results and shows that 
increasing the final irradiance width (diamond, then triangle, then square, then 
X) actually lowered the coherence length in the saturation regime while still 
following the theory in the Rytov regime. This behavior indicated that the 
observed decrease in the saturation regime was physical and not due to 
simulation constraints. 

To further investigate this phenomenon, beam wave (i.e. divergence 
intermediate between spherical and plane waves) and plane wave 
approximations were propagated on a 512x512 grid in which the width of the 
source varied from 4 grid elements (Ax) to 384 grid elements (3/4 of the grid 
width). Figure (64) plots the resulting coherence lengths and indicates that, 
based on the behavior in the Rytov regime, the 4 Ax source (circles) produced 
the spherical wave having large divergence of the E-field, the intermediate 
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Figure 63 Coherence lengths r 0 for varying amounts of spherical divergence, 
hence final irradiance pattern width: diamond = 4/8, triangle = 5/8, square = 6/8, 
X = 7/8 grid width. 
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Figure 64 Coherence lengths for increasing source width on a 512x512 grid 
circles=4Ax (spherical wave); X's=8Ax, squares=16Ax, triangles=32Ax, 
diamonds=64Ax, dots=128Ax, asterisks=256Ax, pluses=384Ax (plane wave). 





source sizes (X=8 Ax, square = 16 Ax, triangle = 32 Ax, diamond = 64 Ax) 
produced beam wave divergences, and wide sources (dot = 128 Ax, asterisk = 
256 Ax, plus = 384 Ax) produced plane waves. 

However, as the turbulence strength approached the saturation regime, 
the behaviors changed. The spherical-type propagation coherence length 
dropped -15% at p 0 2 = 15, the beam wave propagations actually increased in 
coherence length, and the plane wave propagations first decreased -5% before 
increasing -15% at p 0 2 = 15. Some of the unevenness in the beam wave 
coherence lengths occurred because smaller regions of the E-field were used to 
calculated the coherence lengths due to the relatively small size and divergence 
of these waves. 

The cause of these behaviors requires further investigation, but a 
hypothesis can be made. As noted earlier, the spherical wave coherence 
length probably decreased below theory for strong turbulence where the Rytov- 
Tatarski first order perturbation theory was no longer valid. Strong scattering 
may have induced spherical divergence of the beam waves, increasing the 
coherence length, and caused the plane wave approximations to diverge like 
beam waves for high turbulence. 

Recapping the above results, investigations of coherence length via 
computer simulation indicated that first order perturbation theory for coherence 
lengths loses validity in the saturation regime, just as it did for normalized 
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irradiance variance. The behavior of coherence length in the saturation regime 
for spherical, beam, and plane wave cases requires further research. 

These investigations then examined the effect of inner scale upon 
coherence length for a spherically diverging beam, utilizing the 1024x1024 
realizations run for the normalized irradiance variance calculations. Figure (65) 
shows the coherence length r 0 from numerical integration of Eq. (61) for 
Gaussian inner scales of 5, 10, and 15 cm. The Gaussian inner scale made 
the coherence lengths larger by reducing the energy at high spatial frequencies. 
Larger inner scales monotonically produced larger coherence lengths from 
numerical integration at a given turbulence strength (~8% larger at p 0 2 = 1 
(beginning of saturation regime) for Go = 15 cm). A plot of HWHM from 
numerical integration for Gaussian inner scales would appear similar since the 
theoretical HWHM is proportional to the Kolmogorov coherence length r 0 

HWHM = 0.382 r 0 . (147) 

Figure (66) shows the coherence length r 0 and Fig. (67) shows the 
HWHM from wave optics computer simulations with Gaussian inner scales of 5, 
10, and 15 cm (solid lines), superimposed upon the numerical integration 
predictions of Fig. (65) (dotted lines). The r 0 and HWHM simulation values 
agreed well with theory for p 0 2 < 0.1, but started deviating from theory even 
before p 0 2 = 1, i.e. the onset of saturation for the normalized irradiance 
variance. In the saturation regime, computer simulation r 0 and HWHM still 
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Figure 65 Coherence length r 0 from numerical integration for Gaussian inner 
scales of 5, 10, and 15 cm, normalized by the theoretical zero inner scale 
coherence lengths. 




Figure 66 Coherence lengths r 0 from computer simulation with Gaussian inner 
scales of 5, 10, and 15 cm; values normalized by the theoretical zero inner 
scale coherence lengths. 
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Figure 67 HWHM's from computer simulation with Gaussian inner scales of 5 
10, and 15 cm; values normalized by the theoretical zero inner scale HWHM's. 




agreed well with each other, and their behavior showed a combination of 
decrease in coherence length due to saturation regime turbulence and increase 
in coherence length due to inner scale. 

To elucidate the effects of inner scale by itself, Figs. (68) and (69) show 
the same plots of computer simulation coherence length r 0 and HWHM, but now 
normalized by the computer simulation zero inner scale r 0 and HWHM values, 
respectively, effectively removing the saturation contribution. These plots 
clearly show that, even in the saturation regime, inner scale increased the 
coherence length (solid lines) similarly to the predictions of theory (dotted lines). 

The next set of figures plots coherence length and HWHM for 
propagations through turbulence with the Hill viscous-convective enhancement 
inner scale. Figure (70) plots the predicted coherence lengths from numerical 
integration of Eq. (61). Note that the numerical integration coherence lengths 
dropped ~5% in the range p 0 2 = [0.1, 1] before increasing for higher turbulence 
strength. Figures (71) and (72) plot the computer simulation coherence lengths 
r 0 and HWHM, normalized by the theoretical values for spherical waves. In the 
Rytov regime, the coherence length r 0 decreased (~5% for ^ = 15 cm) even for 
very low turbulence (p 0 2 < 0.1 )z, but then increased in the saturation regime as 
for the Gaussian inner scale. Figures (73) and (74) plot the same coherence 
lengths r 0 and HWHM's but divided by the zero inner scale computer simulation 
values to emphasize the effect of inner scale. Figure (73) clearly shows the 
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small decrease in coherence length r 0 in the Rytov regime, and both plots show 
the general agreement between theory and computer simulation for the effect of 
inner scale upon coherence length in the saturation regime. In summary, the 
inner scale increased the coherence length in the saturation regime as much as 
50% compared to the zero inner scale case, and more than compensated for 
the decrease from strong scatter. 
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Figure 68 Coherence lengths r 0 from computer simulation for Gaussian inner 
scales of 5, 10, and 15 cm; values normalized by computer simulation zero 
inner scale r 0 . 
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Figure 69 HWHM's from computer simulation with Gaussian inner scales of 5 
10, and 15 cm; valued normalized by computer simulation zero inner scale 
HWHM's. 
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Figure 70 Coherence length r 0 from numerical integration for Hill viscous- 
convective enhancement inner scales of 5, 10, and 15 cm, normalized by the 
theoretical zero inner scale coherence length. 
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Figure 71 Coherence lengths r 0 from computer simulation with Hill viscous- 
convective enhancement inner scales of 5, 10, and 15 cm; values normalized 
by the theoretical zero inner scale coherence lengths. 
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Figure 72 HWHM's from computer simulation with Hill viscous-convective 
enhancement inner scales of 5, 10, and 15 cm; values normalized by the 
theoretical zero inner scale HWHM's. 
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Figure 73 Coherence lengths r 0 from computer simulation for Hill viscous- 
convective enhancement inner scales of 5, 10, and 15 cm; values normalized 
by computer simulation zero inner scale r 0 . 




HWHM / 



Figure 74 HWHM's from computer simulation with Hill viscous-convective 
enhancement inner scales of 5, 10, and 15 cm; valued normalized by computer 
simulation zero inner scale HWHM's. 





V. CONCLUSIONS 


Variations in the index of refraction within a turbulent medium alter an 
E-field propagating through the medium. The Rytov-Tatarski, perturbation 
theory predicts the effects of the turbulence upon the irradiance statistics and 
coherence length of the propagating E-field. Computer simulations modeled the 
propagation of the E-field through the turbulent medium, producing irradiance 
and coherence statistics to compare with theoretical results. 

These investigations used a split-step, Huygens-Fresnel, wave optics, 
computer simulation to model an E-field propagating through a turbulent 
stratosphere. The limits of validity of the simulations were determined based 
upon aliasing considerations, choice of source, and robustness of statistical 
calculations, and produced the following guidelines: 

• The element size for an NxN grid should satisfy Ax = va L / N . 

• The maximum turbulence strength for which an NxN grid produces valid 
E-fields is given by (3 0 2 max « 0.1 N 09 . 

• A coherence length r 0 * 2.5 Ax corresponds to the maximum turbulence 
strength for which an NxN grid produces valid E-fields. 

• The number of phase screens should be > 30 . 

• The number of realizations should be > 30 . 

• Low spatial frequency corrections to phase screens improve the accuracy 
of the normalized irradiance variance by 5% and of the E-field coherence 
length by 30%. 



• Half width at half maximum of the atmospheric MTF and an iterative fit r 0 
provide the most stable parameterizations of the E-field coherence length. 

• Telltale signs of aliasing include a fine-grained irradiance pattern, a boxed 
perimeter of the irradiance pattern, and peaking of energy toward the 
center of the computation grid. 

This research investigated the effect of (1) zero inner scale, (2) Gaussian 
inner scale, (3) Hill's and (4) Frehlich's viscous-convective enhancement inner 
scales, and (5) grid cutoff inner scale on the normalized irradiance variance of a 
spherical wave propagating through a turbulent medium. For the Rytov regime, 
the normalized irradiance variances with grid cutoff, Gaussian, and Hill/Frehlich 
viscous-convective enhancement inner scales were compared to the zero inner 
scale case and to the values from numerical integration of the Rytov-Tatarski 
predictions. For low turbulence strengths, the variances obtained from the 
simulations agreed within 2% of the values from the numerical integrations. 

The grid cutoff inner scale, implicit in discrete grid wave optics computer 
simulations, affected the variance negligibly compared to a true zero inner scale 
at low turbulence strengths with the large 1024x1024 grid. Application of a 
Gaussian inner scale reduced the normalized irradiance variance as much as 
40% for small p 0 2 compared to the zero inner scale case. The more realistic 
Hill viscous-convective enhancement inner scale raised the normalized 
irradiance variance by up to 30% for smaller inner scale values, but for larger 
values of inner scale eventually reduced the variance below the zero inner 
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scale value as much as 30%. The latter contrasted with the behavior in the 
saturation regime (p 0 2 > 1) where larger inner scales continually enhanced the 
normalized irradiance variance. The Frehlich parameterization of the viscous- 
convective enhancement gave normalized irradiance values that agreed within 
3% of the Hill inner scale values over the entire range of turbulence strengths 
investigated. 

The coherence of the E-field was studied by computing the average 
atmospheric MTF from the propagated fields. Parameterizing the MTF with a 
coherence length r 0 and half width at half max (HWHM) allowed comparison of 
the coherence length of the E-field with the predicted coherence length from the 
Rytov-Tatarski-Fried theory. In the Rytov regime, simulation coherence lengths 
and HWHM's for spherical and plane wave approximations agreed within 5% of 
the theoretical coherence lengths/HWHM's for zero inner scale. However, in 
the saturation regime, the spherical wave coherence length decreased as much 
as 25% below the theory. Similar decreases resulted for different widths of the 
final E-field. Beam wave approximations gave coherence lengths that 
increased toward the spherical wave values in the saturation regime, while 
plane wave approximations deviated from ~5% below to -15% above the theory 
in the saturation regime. Increasing inner scale increased the coherence length 
of a spherically diverging E-field by up to 50% relative to the zero Inner scale 
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case in the saturation regime. The amount of the increase agreed with 
numerical predictions from the analytic theory. 

Several avenues for further research exist. Larger grid sizes and finer 
mesh could explore the behavior of normalized irradiance variance and 
coherence length at higher turbulence strengths. In particular, the behavior of 
coherence length for spherical, beam, and plane waves in strong turbulence 
conditions need a more precise description. The atmospheric MTF's from 
simulation and theory could be compared directly, rather than through a 
coherence length parameterization, and the structure functions could also be 
calculated and compared directly to investigate whether the 2/3 power law 
structure function holds in saturation. Simulations could allow C n 2 (z) to vary 
along the path and could include an outer scale in the spectrum of refractive 
index fluctuations to study the effects of large scale anisotropy of atmospheric 
turbulence on normalized irradiance variance and coherence length (possibly 
basing the C n 2 (z) variations upon high frequency radar data that can resolve the 
large scale variations down to about 300 m). Larger ( > 1 Gbyte RAM) and 
faster computing resources would provide the catalyst for all such further 
investigations of the propagation of an E-field through turbulence. 
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APPENDIX 


Sample listing of YAPS event input file: 


EVENTS: 


EXPLANATION: 


t 7 idebug flag and random number seed 

0.0 0.053033 ;zenith angle and r_0 at 0.5 microns 

2 ;number of wavelengths 

0.5e-6 0.5e-6 ;list of wavelengths 

'surf 'beam' 0 0 0 2.531058298 0.0 0.0 0.00 

;surf name, vertex location, clear aperture, 
;super-Gaussian exponent, inner scale 

'surf 'atmosl' 0 0 6250 8.0 0.0 0.0 0.00 

;surf name, vertex location, clear aperture, 
;super-Gaussian exponent, inner scale 
'prof 1024 1024 0.009886946 'none' 'dummy' 

;surface size, grid element size, file flags 
'end' ;end of surface summary 

2 1024 1024 ;number of fields and dimensions 

'times' 0.0 0.090002 ;time initialization 

'thread' 0.0 1.0 1 .propagation start 

'finit* 1 0.009886946 2 0 0 0 1 0 0 200000 +1 

; initialize field 
; apply aperture profile 
;convert from amplitude/phase to complex 
;change focus to apply spherical phase 
;back propagate to create source 
;copy field to use as source later 
;propagate 
;open field output file 
;remove spherical phase from field 
;save field to output file 
;copy source field to working grid 


'apsrf 1 1 'beam' 

’aptou’ 1 

'chgfcs' 1 0 0 1.0e+30 
'prop 1 1 0 0 200000 
'fldcp' 1 2 
'prop' 2 0 0 0 

'openff '/data/davis/cl 18' 11 
'chgfcs' 2 0 0 200000 
'svflddx' 2 11 385 640 385 640 
'fldcp' 1 2 
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(following three steps repeated 32 times to propagate a distance L ) 

'prop' 2 0 0 193750 ;propagate field 

'mkscrn' 6250.0 1.0e-18 193750.0 ’atmosl' 

;create phase screen for Az, Cn2, position 
'apsrf 2 1 'atmosl' ;apply phase screen to field 

( save the field ) 

'chgfcs' 2 0 0 200000 ;remove spherical phase from field 

'svflddx' 2 11 385 640 385 640 ;save field to output file 

'fldcp' 1 2 '.copy source field to working grid 

( repeat above propagation 30 times ) 


'closefl' 11 ;close output file 

'end' ;end simulation 
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